生成表达矩阵

生信喵

背景
       RNAseq 分析的核心是生成一个表达矩阵，该矩阵结构中，行为基因名或转录本名，列为样品名。
      
       RNAseq 分析流程图

一、质控过滤
1.1 fastqc 质控数据

1. # fastqc
   
2. mkdir qc
   
3. fastqc -t 16 -f fastq -o qc ./fastq/*.gz
   
4. #multiqc整合报告
   
5. multiqc -d qc/ -o multiqc

复制代码

1.2 fastp 过滤

1. ls -1 *.gz | xargs -n 2 | while read {i,j};do echo "fastp -i ${i} -I ${j} -o ../clean/${i%_*}_clean.1.fq.gz -O ../clean/${i%_*}_clean.2.fq.gz -z 4 -q 20 -u 30 -n 10 -h ../clean/${i%_*}.clean.html" >>fastp.sh;done;
   
2. parallel -j 8 -a fastp.sh

复制代码

二、与参考序列比对
       将测序数据与参考序列的比对工具有很多选择，例如 bwa，bowtie2，minimap2，tophat2，hisat2，STAR，subread 等很多工具，由于可变剪切的存在需要对 reads 进行 splice 比对，因此对于 RNAseq 分析需要选择适合 splice比对的工具，这里推荐使用 hisat2 与 STAR，如果是三代测序 RNAseq 数据，可以选择 minimap2 或者 STARlong 工具。

       39 个 RNAseq 比对工具比对的文章：

1. https://www.nature.com/articles/s41467-017-00050-4

复制代码

      比对之后得到的 bam 文件可以使用 featureCounts 或者 HTSeq 工具进行 reads计数。这四种工具可以组合使用，例如比对可以选择 hisat2 与 STAR，计数可以使用 featureCounts 与 HTSeq。
       这里选择两种方案用于比较，当然也可以选择 STAR 与 featureCounts 的方案。

三、hisat2+featureCounts 比对方案
3.1 hisat2 建立索引

1. #方法一：下载索引
   
2. # https://benlangmead.github.io/aws-indexes/hisat
   
3. wget https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz
   
4. 
   
5. #方法二：自建索引
   
6. #演示以前面提取的21号染色体为例
   
7. #1 提取外显子信息
   
8. hisat2_extract_exons.py chr21.gtf >chr21.exon
   
9. 
   
10. #2 提取剪切位点
    
11. hisat2_extract_splice_sites.py chr21.gtf >chr21.ss
    
12. 
    
13. #3 构建索引
    
14. hisat2-build -p 12 --exon chr21.exon --ss chr21.ss chr21.fa chr211.fa
    
15. #chr211.fa即为构建好的索引，速度很快

复制代码

3.2hisat2比对

1. #hisat2比对,以前面构建的chr21为索引
   
2. ls -1 *_clean.*.fq.gz | xargs -n 2 | \
   
3.     while read {i,j};do \
   
4.     echo "hisat2 -p 4 -x /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr211.fa \
   
5.     -1 ${i} -2 ${j} -S ${i%_*}.sam >${i%_*}.log 2>${i%_*}.err; \
   
6.     /share/home/xiehs/tools/samtools/samtools-1.10/samtools sort -@ 4 -o ${i%_*}.sorted.bam ${i%_*}.sam;\
   
7.     /share/home/xiehs/tools/samtools/samtools-1.10/samtools index ${i%_*}.sorted.bam;" >>hisat2.sh\
   
8.     ;done;
   
9.    
   
10. echo "parallel -j 8 -a hisat2.sh" >para.sh
    
11. #集群使用32线程
    
12. #bsub -q fat -n 32 -o %J.log -e %J.err sh para.sh

复制代码

3.3 featureCount 表达矩阵

1. featureCounts -g gene_id -a /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr21.gtf --primary -T 10 -o CountMatrix.txt *.bam
   
2. #去除2-6列
   
3. cut -f 1,7-14 CountMatrix.txt |grep -v "#" >CountMatrix.tsv
   
4. #替换样品名
   
5. sed -i 's#.sorted.bam##g' CountMatrix.tsv

复制代码

四、STAR+HTSeq 比对方案

4.1 STAR 建立索引

1. #方法一：从10x官网下载索引
   
2. #https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build#GRCh38_2020A
   
3. wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
   
4. 
   
5. #方法二：构建索引 STAR
   
6. STAR --runThreadN 12 --runMode genomeGenerate --genomeDir star \
   
7.     --genomeFastaFiles chr21.fa \
   
8.     --sjdbGTFfile chr21.gtf

复制代码

4.2 STAR 比对

1. ls -1 *_clean.*.fq.gz | xargs -n 2 | \
   
2.     while read {i,j};do \
   
3.     echo "STAR --genomeDir /share/home/xiehs/14.rnaseq/ref/ENSEMBL/star --runThreadN 6 \
   
4.     --readFilesIn ${i} ${j} --readFilesCommand zcat \
   
5.     --outFileNamePrefix ${i%_*} --outSAMtype BAM Unsorted \
   
6.     --outSAMattributes All 1>${i%_*}.log 2>${i%_*}.err;" >>star.sh\
   
7.     ;done;
   
8. 
   
9. echo "parallel -j 8 -a star.sh" >para1.sh
   
10. #集群使用48线程
    
11. #bsub -q fat -n 48 -o %J.log -e %J.err sh para1.sh

复制代码

4.3 htseq 计数

1. ls -1 *.bam | while read i;do echo \
   
2.     "htseq-count -f bam -t exon -r pos -s no -c ${i%Aligned.out.bam}.tsv -n 12 \
   
3.     ${i} /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr21.gtf 2>${i%Aligned.out.bam}.log"\
   
4.     ;done;

复制代码

生信喵

因为star占用内存大，我们采用hisat2，可以尝试全染色体的gtf比对再定量。最后得到的就是count矩阵了。