去rRNA可以解决GC含量双峰的右峰

前些天我在生信技能树提出来了一个转录组数据分析的疑难杂症：RNA-seq的fastq文件里面为什么有gc含量的双峰，就是fastq测序数据质量控制的时候发现了GC含量的双峰，然后我简单分析了那些高重复的reads，发现一些是BCA文库，一些是rRNA相关的，所以想到了可以首先去除fastq文件中的rRNA ，就安排学徒去探索一下，过程如下：

1. 从NCBI下载rRNA的fa序列

1.1 打开NCBI，选择“Taxonomy” ，搜索 “Homo“

1.2 点击 Homo

1.3 点击 Homo sapiens

1.4 再次点击Homo sapiens

1.5 右边表格点击Nucleotide的“Subtree links”

1.6 左侧选择rRNA

1.7 下载fasta数据

注意：选择显示200个记录，让所有记录都显示在一页，方便下载全部条目。

另存为：hg38_rRNA.fasta，这个文件很小，如果大家完成上面的步骤有困难，也可以留言自己的邮箱，我们直接发送这个文件给大家哈。

1. bowtie2建立rRNA索引

bowtie2-build hg38_rRNA.fasta hg38_rRNA

这个时候的比对工具的选择，并不一定要bowtie2软件哈。

1. 使用bowtie2去除rRNA，重点--un-conc-gz参数。查阅hisat2的帮助文档，发现有同样的参数，所以可以用hisat2完成同样的操作。

index=/data/reference/index/bowtie2/hg38_rRNA/hg38_rRNA
bowtie2 -x ${index} --un-conc-gz SRR6236728_rmrRNA.fq.gz -1 SRR6236728_1_val_1.fq.gz -2 SRR6236728_2_val_2.fq.gz -p 6 | samtools sort -o tmp.bam -
rm tmp.bam

bowtie2运行耗时20min左右，如果是hisat2，可能会更快，大家可以自行测试和比较一下。

1. 去除前后的fq文件比较

上面fq文件包含val标签的是去除rRNA之前的fq文件，包含rmrRNA标签的是去除rRNA之后的fq文件。

下面是不同fq文件的fastqc报告：

可以看到，少了8百万条reads。

可以看到，GC含量最后一个峰是由rRNA导致，因为8百万条reads被去除后，该峰就消失了。

最后剩下的问题，就是GC含量的另外一个峰。我们后续再谈它！