上一篇文章我们提到作者有非常好的代码资源，但是数据没有权限下载。因此，我又重新找了一篇文献，使用这个文章的数据来进行m6A图文复现。数据相关文献简要介绍如下：
例子来自于发表exomePeak这个软件的文献中提供，相应文献为：doi:10.1038/nn.3449，于30 June 2013发表在nature NEUROSCIENCE上，比较老的数据了。现在m6A的实验方法也早已经更新迭代，得到的测序数据也一般以链特异性、双端150bp为主流。
image-20210114211643359.png
1 数据背景
使用Fto缺陷和野生型的小鼠中脑组织，每个表型三个生物学重复进行MeRIP-Seq。
FTO，也称为ALKBH9，是一个去甲基化酶，属于α-KG依赖的加双氧酶ALKB家族蛋白。FTO最初受到关注，是在GWAS研究中，发现它与肥胖相关【Science. Jun 1；316(5829):1341-5，Science. May 11；316(5826):889-94】，然而它的具体功能，作用底物等，一直是未知的。2011年12月，何川教授研究组在Nature Chemical Biology发表研究N6-Methyladenosine in nuclear rna is a majorsubstrate of the obesity-associated FTO，发现FTO主要定位在细胞核中，是RNA的m6A修饰的去甲基化酶。FTO是第一个被发现的RNA去甲基化酶。
这个数据的文章依然有m6A领域大佬的身影：Samie R Jaffrey ，以及与他同一个学校和单位的Kate D Meyer。他有一篇非常有名的文章大家可以去看看，主要描绘了m6A在mRNA上的分布以及特征：Meyer, K.D. et al. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell 149, 1635–1646 (2012) 。
同一时期开创了m6A研究热潮的还有两个大佬：Dan Dominissini 和 Sharon Moshitch-Moshkovitz，这两个人在2012年，2013年发表了两篇文章，第一次从转录组水平上，大范围、高通量地鉴定了人和小鼠m6A的甲基化水平，这两篇文章是：Nature. 2012 Apr 29;485(7397):201-6 和 doi:10.1038/nprot.2012.148
2 数据下载
得到数据编号GSE47217，ENA数据库使用：PRJNA205149 ，再用ascp下载。
# 从ENA数据库得到fastq下载链接
cat -A fastq.url
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866991/SRR866991.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866992/SRR866992.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866993/SRR866993.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866994/SRR866994.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866995/SRR866995.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866996/SRR866996.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866997/SRR866997.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866998/SRR866998.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR866/SRR866999/SRR866999.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR867/SRR867000/SRR867000.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR867/SRR867001/SRR867001.fastq.gz$
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR867/SRR867002/SRR867002.fastq.gz$
# 使用ascp下载，下载到你所在的当前目录
key_ssh=~/biosoft/miniconda3/envs/rna/etc/asperaweb_id_dsa.openssh
cat fastq.url |while read id
do
ascp -k 1 -QT -l 300m -P33001 -i $key_ssh era-fasp@${id} ./
done
数据下载完之后有一个非常重要的步骤，就是完整性的检验。
md5文件信息同样来自ENA数据库下载的数据信息表格中获取，处理成以下格式。
# 验证数据的完整性, 第一列为md5值，第二列为文件名，中间为两个空格
# md5.txt内容如下：
cat md5.txt
95293cada49264ef20bae77c61056d6c SRR866991.fastq.gz
7c60f24668df87641185d768342b6667 SRR866992.fastq.gz
c3abffe6eb7c3ea61245e78db9772174 SRR866993.fastq.gz
ab863dbf38c2e639dd5823f58e07e9ee SRR866994.fastq.gz
41326299f4e5d7867758dd233459c1e9 SRR866995.fastq.gz
1c6e230ac4806e0dcad6927c59e5ff96 SRR866996.fastq.gz
21dcd0adde84059c9fef056078a6faa8 SRR866997.fastq.gz
e08620b138fdc7d33a8f72d949ff8ad4 SRR866998.fastq.gz
199618b6325b007213d9bb3896181f14 SRR866999.fastq.gz
7b23cb95ffdf90f23e2784f31ad712de SRR867000.fastq.gz
3d14c4c8fcd95e5fa4a2c379d23f3088 SRR867001.fastq.gz
f6d46a57cf66a73d6ea4e1a81e81b06d SRR867002.fastq.gz
下载了一晚上，终于下载好了，检验结果都ok
# 检验结果都ok
md5sum -c md5.txt > check
cat check
SRR866991.fastq.gz: 确定
SRR866992.fastq.gz: 确定
SRR866993.fastq.gz: 确定
SRR866994.fastq.gz: 确定
SRR866995.fastq.gz: 确定
SRR866996.fastq.gz: 确定
SRR866997.fastq.gz: 确定
SRR866998.fastq.gz: 确定
SRR866999.fastq.gz: 确定
SRR867000.fastq.gz: 确定
SRR867001.fastq.gz: 确定
SRR867002.fastq.gz: 确定
此外：早期的m6A数据测序片段偏短并且大多数都是单端测序。
3 数据质控和过滤
在确保数据完整性之后，我们对原始数据进行一下简单的质量评估。这里主要还是使用fastqc软件进行简单的评估。
mkdir qc
fastqc -t 20 -o qc/ SRR*.fastq.gz
# 使用MultiQc整合FastQC结果
multiqc *.zip
qc结果：
image-20210713104734504.png
查看所有样本得整合报告：multiqc_report.html
数据有一些N：
image-20210713105006818.png
有一个样本接头含量比较高：
image-20210713105040033.png
有两个两本的GC含量稍微有点异常：
image-20210713105115926.png
有两个样本的数据重复率偏高：
image-20210713105224286.png
这就是以上数据的一个简单评估结果，数据整体Q30挺好，就是还有些接头，N碱基含量，GC分布异常等问题。GC异常的问题，我们后面专门再说，数据质量的好坏与前期样本质量，实验环节息息相关。但无论数据是好是坏，我们生信端能做的就是尽量将异常给去掉来保证后面分析的结果。当然，有钱的可以选择重新提取样本进行测序。
然后使用trim_galore进行过滤，数据读长比较短，我们这里保留大于15bp的reads。
mkdir cleandata
# 过滤
ls ../fastq/SRR*gz |while read id
do
trim_galore --phred33 -q 25 -e 0.1 --length 15 --stringency 3 --fastqc --max_n 3 -o ./ $id >${id}.log
done
# 过滤后的质控
multiqc *.zip
数据到这里就得到了一个cleandata，后面就开始进行比对部分分析了。
image-20210713105537752.png
等待更新~
网友评论