在上一期中我們得到了cleandata后，先使用bowtie2與NCBI的rRNA的序列進(jìn)行比對，進(jìn)行了去除rRNA序列的步驟，得到了去除rRNA之后的數(shù)據(jù)如下：

image

一、rRNA比對結(jié)果

我們先來看下rRNA比對的結(jié)果，以此來查看數(shù)據(jù)中rRNA的含量，并且看看之前有兩個(gè)GC含量分布異常的樣本是不是真的rRNA序列含量比較多，GC異常的兩個(gè)樣本（額，第二期教程里面展示的QC結(jié)果里其實(shí)是數(shù)據(jù)過濾之后的QC結(jié)果，數(shù)據(jù)過濾是進(jìn)行接頭、低質(zhì)量的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結(jié)果中GC含量特征分布圖結(jié)果一致。

rRNA數(shù)據(jù)過高，對m6A Peak Calling分析的主要影響是有效數(shù)據(jù)量變少。

二、去除完rRNA之后，接下來就是與參考基因組的比對了

1.參考基因組下載

這里我們使用ensembl數(shù)據(jù)庫的參考基因，下載方式如下

1.進(jìn)入網(wǎng)址：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

如果網(wǎng)速不好，可能的折騰挺久的。

2.建索引

由于m6A數(shù)據(jù)是RNA測序，因此我們這里使用適用于RNA比對的軟件Hisat2進(jìn)行比對，那么，建立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腳本提交后臺運(yùn)行
# index.sh的內(nèi)容為：
hisat2-build Mus_musculus.GRCm39.dna.primary_assembly.fa.gz   Mus_musculus.GRCm39.dna.primary_assembly

# 提交后臺運(yùn)行，注意養(yǎng)成一個(gè)好習(xí)慣，保存運(yùn)行日志方面查看進(jìn)度
nohup sh index.sh >index.sh.log &

索引建立完了之后會生成以下幾個(gè)以ht2l結(jié)尾（大的參考基因組生成的索引）的文件，如果是其他物種比如人后綴會是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.比對

參考基因組已經(jīng)準(zhǔn)備好了，去除rRNA的fq文件也準(zhǔn)備好了，萬事具備，開干！

# 首先處理得到sampleID，也可以直接從ENA的那個(gè)表格中提取出來
ls -1 bowtie2/*gz | cut -d'/' -f 2 |cut -d'.' -f 1 >sampleId.txt

# sampleId.txt內(nèi)容如下，其實(shí)就是在后面的代碼中起到一個(gè)變量的作用
SRR866991
SRR866992
SRR866993
SRR866994
SRR866995
SRR866996
SRR866997
SRR866998
SRR866999
SRR867000
SRR867001
SRR867002

#將以下內(nèi)容寫入到Hisat.sh腳本中
#注意等號兩邊不能有空格，index內(nèi)容要換成自己的正確路徑，index再次強(qiáng)調(diào)為索引前綴
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

# 提交后臺運(yùn)行
nohup sh Hisat.sh >Hisat.sh.log &

運(yùn)行要幾個(gè)小時(shí)，等一晚上吧~

我們下期再更新，下期就是查看比對結(jié)果，以及bam轉(zhuǎn)bw然后進(jìn)行本教程的第一張圖熱圖繪制。