从开始学习转录组分析,都是按照网上的教程别人的数据一步一步的跑,只要是严格按照教程上来,基本都不会出现错误。在整个过程中,也没有想太多,每一步详细来说,是干什么的,每个参数是怎么设的,都没有经过仔细的考虑。
最近拿到了自己的数据,才发现自己之前做的太简单,缺乏深入的理解。以致遇到很多问题,都难以解答。菜鸟的路还很漫长。
这次分享的这个教训就是关于链特异性建库RNA-seq数据,我按照hisat2 -> samtools -> htseq-count的流程,得到count文件在R上用DEseq2进行下游分析。
得到count文件后,我先在Excel中打开一个文件查看,计算了下总得read数(如下图),发现比我比对上的read数少了一半,感觉不太对,却还不知道问题出在哪
然后我用Subread -> featureCounts -> DESeq2这个流程跑了一遍,得到count文件,发现确实少了差不多一半
然后想起是不是链特异性这个参数要特别设置,然后我开始从头开始看,在hisat2中发现到问题所在
后面的htseq-count参数也要改:(--stranded=reverse)
都修改之后:(得到对应到基因上的count数差不多是之前的2倍)
附上我的跑的代码:(有什么错误,还请指点出来)
、、、
for i in `seq 1 10`
do
hisat2 -p 20 --rna-strandness=FR --dta -x ./hisat2_index/rice_tran -1 **-${i}_paired_1.fq.gz -2 **-${i}_paired_2.fq.gz -S **-${i}_paired.sam 2>**-${i}_paired.log
samtools view -@ 20 -S **-${i}_paired.sam -b > **-${i}_paired.bam
samtools sort -@ 20 -n **-${i}_paired.bam -o **-${i}_paired_sorted.bam
htseq-count -s reverse -r pos -f bam --type=exon --idattr=gene_id **-${i}_paired_sorted.bam all.gtf3 >**-${i}_count.txt 2>**-${i}_count.log
done
、、、
菜鸟之路继续前行。。。