? ? 從開始學習轉錄組分析，都是按照網上的教程別人的數據一步一步的跑，只要是嚴格按照教程上來，基本都不會出現錯誤。在整個過程中，也沒有想太多，每一步詳細來說，是干什么的，每個參數是怎么設的，都沒有經過仔細的考慮。

? ?最近拿到了自己的數據，才發(fā)現自己之前做的太簡單，缺乏深入的理解。以致遇到很多問題，都難以解答。菜鳥的路還很漫長。

? ?這次分享的這個教訓就是關于鏈特異性建庫RNA-seq數據，我按照hisat2?-> samtools -> htseq-count的流程，得到count文件在R上用DEseq2進行下游分析。

? 得到count文件后，我先在Excel中打開一個文件查看，計算了下總得read數（如下圖），發(fā)現比我比對上的read數少了一半，感覺不太對，卻還不知道問題出在哪

count文件（htseq-count）

? 然后我用Subread -> featureCounts -> DESeq2這個流程跑了一遍，得到count文件，發(fā)現確實少了差不多一半

count文件（featurecount）

然后想起是不是鏈特異性這個參數要特別設置，然后我開始從頭開始看，在hisat2中發(fā)現到問題所在

hisat2指南中截取

后面的htseq-count參數也要改：（--stranded=reverse）

htseq-count

都修改之后：（得到對應到基因上的count數差不多是之前的2倍）

修改參數后得到的count文件

附上我的跑的代碼：（有什么錯誤，還請指點出來）

、、、

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

、、、

菜鳥之路繼續(xù)前行。。。