m6A图文复现02-数据下载和质控
前面我们分享了 跟着Nature Medicine学MeDIP-seq数据分析，数据和代码都是公开，这个2G的压缩包文件，足以学习3个月，写60篇教程。同时也分享了 全套MeRIP-seq文章图表复现代码，其实MeRIP-seq其实就是RNA水平的又叫做m6a测序。
但很多粉丝留言表示这些英文教程看不懂，数据也很分散，没有中文解说实在是很难跟下来，希望我们出一个手把手系列教程。
这个全套 MeRIP-seq 图表复现代码在GitHub：https://github.com/al-mcintyre/merip_reanalysis_scripts 这个也是接近2G的压缩包！
其实很早以前我就在《生信技能树》发布过教程：新的ngs流程该如何学习（以CUT&Tag 数据处理为例子），提到了我自己是不太可能去把所有的ngs流程全部录制视频的，只能说是更好的传达学习方法给到大家。其实如果你看过我表观组学，比如《ChIP-seq数据分析》 和 《ATAC-seq数据分析》 就会发现其实这个m6A数据处理大同小异的，当然了，肯定是会有一些细微差异是需要注意的。
虽然我没有时间，但是我们的两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，还是有幸招募到了愿意花时间给大家做整理的小伙伴。所以有了这个MeRIP-seq 图表复现交流群。
上一篇文章我们提到作者有非常好的代码资源，但是数据没有权限下载。因此，我又重新找了一篇文献，使用这个文章的数据来进行m6A图文复现。数据相关文献简要介绍如下：
例子来自于发表exomePeak这个软件的文献中提供，相应文献为：doi:10.1038/nn.3449，于30 June 2013发表在nature NEUROSCIENCE上，比较老的数据了。现在m6A的实验方法也早已经更新迭代，得到的测序数据也一般以链特异性、双端150bp为主流。
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 数据下载
如果你的sratoolkit有问题，或者prefetch命令下载sra文件速度太慢，可以参考：使用ebi数据库直接下载fastq测序数据 , 需要自行配置好，然后去EBI里面搜索到的 fq.txt 路径文件。
得到数据编号GSE47217，ENA数据库使用：PRJNA205149 ：
- https://www.ncbi.nlm.nih.gov/bioproject/PRJNA205149
再用ascp下载，其中aspera是提供conda安装的哦！
# 从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
查看所有样本的整合报告：multiqc_report.html
数据有一些N：
image-20210713105006818
有一个样本接头含量比较高：
有一个样本接头含量比较高
有两个样本的GC含量稍微有点异常：
两个样本的GC含量稍微有点异常
有两个样本的数据重复率偏高：
数据质量控制表格
这就是以上数据的一个简单评估结果，数据整体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，后面就开始进行比对部分分析了。
cleandata
等待更新~