m6A图文复现02-数据下载和质控 / 开普饭

但很多粉丝留言表示这些英文教程看不懂，数据也很分散，没有中文解说实在是很难跟下来，希望我们出一个手把手系列教程。

这个全套 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：

有一个样本接头含量比较高：

有两个样本的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，后面就开始进行比对部分分析了。

等待更新~

文末友情推荐

如果你也想进群交流，可以看：m6A-seq 数据处理及图表复现交流群

m6A图文复现02-数据下载和质控

1 数据背景

2 数据下载

3 数据质控和过滤

文末友情推荐

相关推荐