比对

这一步将使用hisat2或STAR对数据进行比对。hisat2在比对速度与资源占用上说比STAR要好上不少，但是也有不少人建议使用STAR，因为STAR可以获得更好的比对结果。因此具体的选择还是根据实际情况定。关于hisat2、STAR、topHat的对比可以看这篇文章。

HISAT2

分步进行

原则上，hisat2先将数据比对到sam文件，再使用samtools转换为bam文件后再进行排序。其中-x后跟的是索引文件的名字（不包括扩展名）。hisat2对SRR1374921比对用时3分钟，比对率在98.24%。

hisat2 -p 10 \
	-x genome/mm39 \
	-U rRNA/SRR1374921_filtered.fq.gz
	-S bam/SRR1374921.sam

再使用samtools转为bam

samtools view -bS bam/SRR1374921.sam \
	| samtools sort -@ 10 - -o bam/SRR1374921.bam

最后再用samtools进行排序并创建索引

samtools index bam/SRR1374921.bam

然后其实可把SRR1374921.sam删掉了。

合并进行

其实我们也可以一条命令，直接生成排序后的bam。使用管道和tee同时生成比对报告。

hisat2 -p 10 \
	-x genome/mm39 \
	-U rRNA/SRR1374921_filtered.fq.gz \
	| tee >(samtools flagstat - > bam/SRR1374921.hisat2.flagstat.txt) \
	| samtools sort -O BAM \
	| tee bam/SRR1374921.bam \
	| samtools index - bam/SRR1374921.bam.bai

写一个循环来跑这些数据。下载的内容里有个mergelist.txt，有所有数据样本名，把样本名截下来。

cut mergelist.txt -d'_' -f1 > samplelist.txt

然后循环运行

cat samplelist.txt | while read line; do \
	hisat2 -p 10 -x genome/mm39 \
	-U rRNA/${line}_filtered.fq.gz \
	| tee >(samtools flagstat - > bam/${line}.hisat2.flagstat.txt) \
	| samtools sort -O BAM \
	| tee bam/${line}.bam \
	| samtools index - bam/${line}.bam.bai; done

比对完成后，进入下一步统计reads数。

STAR

同样的数据，STAR需要50min才能完成，比对率85.76%。

STAR --genomeDir genome/star_index/ \
	--readFilesIn rRNA/SRR1374921_filtered.fq.gz \
	--runThreadN 10 --outSAMtype BAM SortedByCoordinate \
	--quantMode GeneCounts \
	--outFileNamePrefix bam/SRR1374921. \
	--readFilesCommand zcat
mv bam/SRR1374921.Aligned.sortedByCoord.out.bam bam/SRR1374921.star.bam
samtools index bam/SRR1374921.star.bam

对文件存放结构不满意的需要把log移动到log文件夹。

由于在索引中还加入了GTF，STAR会同时输出比对到各个基因的counts，结果在ReadsPerGene.out.tab中，后三列分别是“counts for unstranded RNA-seq”， “counts for the 1st read strand aligned with RNA (htseq-count option -s yes)”， “counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)”。根据需求选用哪一列。