m6A图文复现03-测序数据去除rRNA序列并且比对到参考基因组
m6A图文复现03-测序数据去除rRNA序列并且比对到参考基因组
前面我们分享了 跟着Nature Medicine学MeDIP-seq数据分析，数据和代码都是公开，这个2G的压缩包文件，足以学习3个月，写60篇教程。同时也分享了 全套MeRIP-seq文章图表复现代码，其实MeRIP-seq其实就是RNA水平的，又叫做m6a测序。
下面是MeRIP-seq 图表复现笔记
在上一期：m6A图文复现02-数据下载和质控
中我们得到了cleandata，接下来是要比对到参考基因上进行比对过程的分析。一般来说，在比对之前，我们可以选择先去除rRNA序列然后再与参考基因组进行比对，文章中采用的策略也是如此：
那么你或许会有以下几个随着而来的问题：
- 1.为什么测序数据中会出现rRNA？
- 2.rRNA不去除与去除会对我后面的分析造成什么样的影响？
rRNA背景知识
在进行具体操作之前，我们可以来先了解一下rRNA这个东西。
首先，看rRNA在整个细胞中的组成比例：从mass指标来看，一个哺乳动物的细胞中rRNA占比最多，能达到80-90%
比例：
Note：Estimate of RNA levels in a typical mammalian cell。ref：Front Genet. 2015 Jan 26;6:2
rRNA在参考基因组上存在多个copy（https://www.britannica.com/science/nucleic-acid/Ribosomal-RNA-rRNA#ref594026），比对时会造成数据比对结果中多比率偏高。rRNA的存在对后续的影响就主要看后面的分析使用的是唯一比对reads还是所有比对上的reads了。
那么，数据中为什么会测到rRNA呢？这就跟前期建库采用的策略有关系了，去除rRNA的试剂盒的效率等问题。
此次，我们先去文中提到的网址下载rRNA序列，然后使用bowtie2进行比对。
rRNA序列下载
打开NCBI网址：https://www.ncbi.nlm.nih.gov/，选择Taxonomy数据库，输入10090，10090为小鼠的物种ID
选择右侧 Nucleotide
保存
下载下来的数据：fasta格式
确认一下下载是否完整，为82条：
rRNA比对保存没有比对上的序列
首先，建索引：
bowtie2-build Mus_musculus.rRNA.fa Mus_musculus.rRNA
索引结果：
比对：
mkdir bowtie2
# 换成自己的路径,注意index为索引前缀
index=/path/rRNA_index/Mus_musculus.rRNA
od=/path/bowtie2
# 单端数据
ls *gz |while read id
do
sam=${id%_*};
nohup bowtie2 -x $index --un-gz ${od}/${sam}.derRNA.fq.gz -U $id -p 8 -S ${od}/${sam}.rRNA.mapped.sam 2>${od}/${sam}.log &
done
# 双端数据
可以看文献中的代码，文献中提供的为双端数据
比对完之后的数据：
后续更新~
在上一期中我们得到了cleandata后，先使用bowtie2与NCBI的rRNA的序列进行比对，进行了去除rRNA序列的步骤，得到了去除rRNA之后的数据如下：
一、rRNA比对结果
我们先来看下rRNA比对的结果，以此来查看数据中rRNA的含量，并且看看之前有两个GC含量分布异常的样本是不是真的rRNA序列含量比较多，GC异常的两个样本（额，第二期教程里面展示的QC结果里其实是数据过滤之后的QC结果，数据过滤是进行接头、低质量的reads过滤，是不能过滤掉rRNA序列的）：
rRNA比对率如下：
cd bowtie2
ls SRR*log >ID1
cat SRR*log |grep 'overall alignment rate' >ID2
paste -d'\t' ID1 ID2
SRR866991.log 25.24% overall alignment rate
SRR866992.log 70.00% overall alignment rate
SRR866993.log 28.12% overall alignment rate
SRR866994.log 91.21% overall alignment rate
SRR866995.log 37.59% overall alignment rate
SRR866996.log 78.50% overall alignment rate
SRR866997.log 35.99% overall alignment rate
SRR866998.log 89.05% overall alignment rate
SRR866999.log 34.60% overall alignment rate
SRR867000.log 62.99% overall alignment rate
SRR867001.log 34.18% overall alignment rate
SRR867002.log 68.39% overall alignment rate
我们可以看到SRR866994和SRR866998的总比对率最高，分别为91.21%和89.05%，与我们观察到的qc结果中GC含量特征分布图结果一致。
rRNA数据过高，对m6A Peak Calling分析的主要影响是有效数据量变少。
二、去除完rRNA之后，接下来就是与参考基因组的比对了
1.参考基因组下载
这里我们使用ensembl数据库的参考基因，下载方式如下
1.进入网址：http://ftp.ensembl.org/pub/release-104/fasta/mus_musculus/dna/
选择文件：Mus_musculus.GRCm39.dna.primary_assembly.fa.gz
下载命令：
mkdir -p Genome/Mouse/GRCm38.104
wget -c http://ftp.ensembl.org/pub/release-104/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz
如果网速不好，可能的折腾挺久的。
2.建索引
由于m6A数据是RNA测序，因此我们这里使用适用于RNA比对的软件Hisat2进行比对，那么，建立Hisat2的参考基因组索引为：
# 注意Homo_sapiens.GRCh38.dna.primary_assembly为索引前缀
hisat2-build Mus_musculus.GRCm39.dna.primary_assembly.fa.gz \ Mus_musculus.GRCm39.dna.primary_assembly
# 建立索引的过程大约为30分钟左右，可以写入sh脚本提交后台运行
# index.sh的内容为：
hisat2-build Mus_musculus.GRCm39.dna.primary_assembly.fa.gz Mus_musculus.GRCm39.dna.primary_assembly
# 提交后台运行，注意养成一个好习惯，保存运行日志方面查看进度
nohup sh index.sh >index.sh.log &
索引建立完了之后会生成以下几个以ht2l结尾（大的参考基因组生成的索引）的文件，如果是其他物种比如人后缀会是ht2：
Mus_musculus.GRCm39.dna.primary_assembly.1.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.2.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.3.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.4.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.5.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.6.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.7.ht2l
Mus_musculus.GRCm39.dna.primary_assembly.8.ht2l
3.比对
参考基因组已经准备好了，去除rRNA的fq文件也准备好了，万事具备，开干！
# 首先处理得到sampleID，也可以直接从ENA的那个表格中提取出来
ls -1 bowtie2/*gz | cut -d'/' -f 2 |cut -d'.' -f 1 >sampleId.txt
# sampleId.txt内容如下，其实就是在后面的代码中起到一个变量的作用
SRR866991
SRR866992
SRR866993
SRR866994
SRR866995
SRR866996
SRR866997
SRR866998
SRR866999
SRR867000
SRR867001
SRR867002
#将以下内容写入到Hisat.sh脚本中
#注意等号两边不能有空格，index内容要换成自己的正确路径，index再次强调为索引前缀
index=Genome/Mouse/GRCm38.104/Mus_musculus.GRCm39.dna.primary_assembly
indir=./bowtie2
outdir=./mapping
cat sampleId.txt | while read id
do
hisat2 -p 10 -x ${index} -U ${indir}/${id}.derRNA.fq.gz 2>${outdir}/${id}.log | samtools sort -@ 10 -o ${outdir}/${id}.Hisat_aln.sorted.bam - && samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai
done
# 提交后台运行
nohup sh Hisat.sh >Hisat.sh.log &
运行要几个小时，等一晚上吧~
我们下期再更新，下期就是查看比对结果，以及bam转bw然后进行本教程的第一张图热图绘制。