生成表达矩阵
本帖最后由 生信喵 于 2022-10-9 17:42 编辑背景
RNAseq 分析的核心是生成一个表达矩阵，该矩阵结构中，行为基因名或转录本名，列为样品名。
RNAseq 分析流程图
一、质控过滤
1.1 fastqc 质控数据
# fastqc
mkdir qc
fastqc -t 16 -f fastq -o qc ./fastq/*.gz
#multiqc整合报告
multiqc -d qc/ -o multiqc1.2 fastp 过滤
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;
parallel -j 8 -a fastp.sh
二、与参考序列比对
将测序数据与参考序列的比对工具有很多选择，例如 bwa，bowtie2，minimap2，tophat2，hisat2，STAR，subread 等很多工具，由于可变剪切的存在需要对 reads 进行 splice 比对，因此对于 RNAseq 分析需要选择适合 splice比对的工具，这里推荐使用 hisat2 与 STAR，如果是三代测序 RNAseq 数据，可以选择 minimap2 或者 STARlong 工具。
39 个 RNAseq 比对工具比对的文章：
https://www.nature.com/articles/s41467-017-00050-4 比对之后得到的 bam 文件可以使用 featureCounts 或者 HTSeq 工具进行 reads计数。这四种工具可以组合使用，例如比对可以选择 hisat2 与 STAR，计数可以使用 featureCounts 与 HTSeq。
这里选择两种方案用于比较，当然也可以选择 STAR 与 featureCounts 的方案。
三、hisat2+featureCounts 比对方案
3.1 hisat2 建立索引
#方法一：下载索引
# https://benlangmead.github.io/aws-indexes/hisat
wget https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz
#方法二：自建索引
#演示以前面提取的21号染色体为例
#1 提取外显子信息
hisat2_extract_exons.py chr21.gtf >chr21.exon
#2 提取剪切位点
hisat2_extract_splice_sites.py chr21.gtf >chr21.ss
#3 构建索引
hisat2-build -p 12 --exon chr21.exon --ss chr21.ss chr21.fa chr211.fa
#chr211.fa即为构建好的索引，速度很快
3.2hisat2比对
#hisat2比对,以前面构建的chr21为索引
ls -1 *_clean.*.fq.gz | xargs -n 2 | \
while read {i,j};do \
echo "hisat2 -p 4 -x /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr211.fa \
-1 ${i} -2 ${j} -S ${i%_*}.sam >${i%_*}.log 2>${i%_*}.err; \
/share/home/xiehs/tools/samtools/samtools-1.10/samtools sort -@ 4 -o ${i%_*}.sorted.bam ${i%_*}.sam;\
/share/home/xiehs/tools/samtools/samtools-1.10/samtools index ${i%_*}.sorted.bam;" >>hisat2.sh\
;done;
echo "parallel -j 8 -a hisat2.sh" >para.sh
#集群使用32线程
#bsub -q fat -n 32 -o %J.log -e %J.err sh para.sh
3.3 featureCount 表达矩阵
featureCounts -g gene_id -a /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr21.gtf --primary -T 10 -o CountMatrix.txt *.bam
#去除2-6列
cut -f 1,7-14 CountMatrix.txt |grep -v "#" >CountMatrix.tsv
#替换样品名
sed -i 's#.sorted.bam##g' CountMatrix.tsv
四、STAR+HTSeq 比对方案
4.1 STAR 建立索引
#方法一：从10x官网下载索引
#https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build#GRCh38_2020A
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
#方法二：构建索引 STAR
STAR --runThreadN 12 --runMode genomeGenerate --genomeDir star \
--genomeFastaFiles chr21.fa \
--sjdbGTFfile chr21.gtf
4.2 STAR 比对
ls -1 *_clean.*.fq.gz | xargs -n 2 | \
while read {i,j};do \
echo "STAR --genomeDir /share/home/xiehs/14.rnaseq/ref/ENSEMBL/star --runThreadN 6 \
--readFilesIn ${i} ${j} --readFilesCommand zcat \
--outFileNamePrefix ${i%_*} --outSAMtype BAM Unsorted \
--outSAMattributes All 1>${i%_*}.log 2>${i%_*}.err;" >>star.sh\
;done;
echo "parallel -j 8 -a star.sh" >para1.sh
#集群使用48线程
#bsub -q fat -n 48 -o %J.log -e %J.err sh para1.sh
4.3 htseq 计数
ls -1 *.bam | while read i;do echo \
"htseq-count -f bam -t exon -r pos -s no -c ${i%Aligned.out.bam}.tsv -n 12 \
${i} /share/home/xiehs/14.rnaseq/ref/ENSEMBL/hisat2/chr21.gtf 2>${i%Aligned.out.bam}.log"\
;done; 因为star占用内存大，我们采用hisat2，可以尝试全染色体的gtf比对再定量。最后得到的就是count矩阵了。
页:
[1]