Skip to content

比对

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

HISAT2

分步进行

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

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

再使用samtools转为bam 

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

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

1
samtools index bam/SRR1374921.bam

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

合并进行

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

1
2
3
4
5
6
7
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，有所有数据样本名，把样本名截下来。 

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

然后循环运行 

1
2
3
4
5
6
7
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%。

1
2
3
4
5
6
7
8
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)”。根据需求选用哪一列。