对于RNA-seq实验与分析流程是三天前开始学习的,简单记录一下。
RNA-seq 实验可以捕获全部RNA(不区分类型),也可以根据成熟转录本尾部特异性polyA 尾巴特异性选择mRNA。
普通的RNA-seq是不能区分链的,也就是说我们不能知道转录本来自正链还是负链,但是通过dUTP的掺入,可用特定的酶将反转录合成的第二条cDNA链降解,这样就知道转录本来自于哪一条链,后续比对到参考基因的时要用特定的参数rna-strandness。
数据处理的过程主要包括质控——映射到参考基因组——组装定量——后续分析(如:差异表达)
1.质控我是用fastp 主要去接头污染和低质量的片段,碱基。
2.质控后用fastqc看一眼,各项指标是否满足向下分析的要求
3.mapping 我使用的是hisat2
(1)用hisat2要提前准备好index,自己建立index比较费时我选择在hisat2官网直接下载
Download | HISAT2 (daehwankimlab.github.io)
(2)index建立过程是需要gtf文件的,这个要留意hisat2官网使用的gtf版本,因为后续组装的时候需要用到gtf文件,要使两者保持一致,如果后续使用ballgown 要加上-dta参数
GENCODE - Human Release 19 (gencodegenes.org)
4.对数据进行定量
(1)stringtie
stringtie -G ../reference/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf -o SRR{i}.gtf -l SRR${i} SRR${i}.q10.sort.bam
ls -l SRR*.gtf | awk '{print $9}' > sample_assembly_gtf_list.txt
stringtie --merge -p 8 -o stringtie_merged.gtf -G ../reference/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf sample_assembly_gtf_list.txt
gffcompare -r ../reference/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf -o gffcompare stringtie_merged.gtf
(2)featureCounts
featureCounts -p -t exon -g gene_id -M -T 8 -a ../reference/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf -o all.featurecounts.txt ./*sort.bam
awk -F '\t' '{print $1,$7,$8,$9,$10,$11,$12}' OFS='\t' all.featurecounts.txt > all_fcount.matrix.txt
5.差异表达分析
使用R包ballgown和DEseq2都可以