m6ASeq RNA甲基化测序

m6A图文复现04-rRNA比对结果与Hisat2比对

2022-01-06  本文已影响0人  信你个鬼

在上一期中我们得到了cleandata后,先使用bowtie2与NCBI的rRNA的序列进行比对,进行了去除rRNA序列的步骤,得到了去除rRNA之后的数据如下:

image

一、rRNA比对结果

我们先来看下rRNA比对的结果,以此来查看数据中rRNA的含量,并且看看之前有两个GC含量分布异常的样本是不是真的rRNA序列含量比较多,GC异常的两个样本(额,第二期教程里面展示的QC结果里其实是数据过滤之后的QC结果,数据过滤是进行接头、低质量的reads过滤,是不能过滤掉rRNA序列的):

image

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然后进行本教程的第一张图热图绘制。

上一篇下一篇

猜你喜欢

热点阅读