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

上一篇文章我们提到作者有非常好的代码资源，但是数据没有权限下载。因此，我又重新找了一篇文献，使用这个文章的数据来进行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

等待更新~