MeRIP-seq 数据分析之数据下载

点击上方关注“公众号”

1背景

关于什么是 m6A 和 MeRIP-seq 这里我就不说了,具体参考：Introduction of m6A。今天我们来分析最近一篇文章的 MeRIP-seq 数据。这篇文章 2021 年 9 月 30 日发表在 nature methods 上：

这篇文章的数据上传到了 GEO 数据库，编号为 GSE151028，样本有 98 个！我们下载其中几个练习就行了，样本如下：

• 总共 10 个，分别为 IVT 的 IP 和 Input，METTL3-KO 的 IP 和 Input 及野生型 WT 的 IP 和 Input，分别有 2 个生物学重复，IVT 就下载一个生物学重复就行了。
• 物种为小鼠。
• 文库类型为双端 150bp 测序，链特异性测序（dUTP 法）。

---

文章主要讲的是：

通过体外合成的没有修饰的 RNA（m6A、m5C 等）来做同样的 MeRIP-seq 测序（称为IVT-seq），相比于正常的 MeRIP-seq 测序，在 IVT-seq 中鉴定的 m6A peak 作为一个假阳性的结果，把正常的 MeRIP-seq 测序的 m6A peak 去除掉这些假阳性结果可以大大提高 m6A 鉴定的准确性。

何川教授 给予点评：这个方法可能将来称为 MeRIP-seq 测序的金标准！

下面是 体外合成没有修饰的 RNA 原理示意图（体外转录实验）：

然后是过滤假阳性 peak 的示意图：

2前期准备

高通量数据下载方式很多，可以使用 prefetch 指定 SRR 号下载 sra 文件、可以去 ENA 数据库下载 fastq 文件，也可以使用 aspera 高速下载、迅雷 也可以，kingfisher 可以指定从多个来源下载，包括以上的。

这里尝试一下 kingfisher 来下载。

安装：

# 下载
$ git clone https://github.com/wwood/kingfisher-download.git
# 进入文件夹
$ cd kingfisher-download/
# 创建环境并安装依赖软件
$ conda env create -n kingfisher -f kingfisher.yml

# 激活环境
$ conda activate kingfisher
# 添加到环境变量
$ cd bin
$ export PATH=$PWD:$PATH

# 查看帮助文档
$ kingfisher -h

                ...::: Kingfisher v0.0.1-dev :::...

  get         -> Download and extract sequence data from SRA or ENA
  extract     -> extract .sra format files
  annotate    -> annotate runs by their metadata

  Use kingfisher <command> -h for command-specific help.
  Some commands also have an extended --full-help flag.

下载数据主要有以下几个来源：

使用第一种 ena-ascp 是最快的，但可能会报错：

$ kingfisher get -r ERR1739691 -m ena-ascp
10/10/2021 08:59:33 AM INFO: Kingfisher v0.0.1-dev
10/10/2021 08:59:33 AM INFO: Attempting download method ena-ascp ..
10/10/2021 08:59:33 AM INFO: Using aspera ssh key file: $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh
10/10/2021 08:59:33 AM INFO: Querying ENA for FTP paths for ERR1739691..
10/10/2021 08:59:35 AM INFO: Downloading 2 FTP read set(s): ftp.sra.ebi.ac.uk/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_1.fastq.gz, ftp.sra.ebi.ac.uk/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_2.fastq.gz
10/10/2021 08:59:35 AM INFO: Running command: ascp -T -l 300m -P33001 -k 2 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh [e-mail]:/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_1.fastq.gz .
10/10/2021 08:59:35 AM WARNING: Error downloading from ENA with ASCP: Command ascp -T -l 300m -P33001 -k 2 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh [e-mail]:/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_1.fastq.gz . returned non-zero exit status 1.
STDERR was: b'ascp: Private key file not found at path /home/junjun/.aspera/connect/etc/asperaweb_id_dsa.openssh, exiting.\n'STDOUT was: b'\r\nSession Stop  (Error: Private key file not found at path /home/junjun/.aspera/connect/etc/asperaweb_id_dsa.openssh)\n'
10/10/2021 08:59:35 AM WARNING: Method ena-ascp failed
Traceback (most recent call last):
  File "/home/junjun/biosoft/kingfisher-download/bin/kingfisher", line 280, in <module>
    main()
  File "/home/junjun/biosoft/kingfisher-download/bin/kingfisher", line 232, in main
    kingfisher.download_and_extract(
  File "/home/junjun/biosoft/kingfisher-download/bin/../kingfisher/__init__.py", line 37, in download_and_extract
    download_and_extract_one_run(run, **kwargs)
  File "/home/junjun/biosoft/kingfisher-download/bin/../kingfisher/__init__.py", line 276, in download_and_extract_one_run
    raise Exception("No more specified download methods, cannot continue")
Exception: No more specified download methods, cannot continue

解决办法：

# 安装ascp
$ conda install -c hcc aspera-cli -y
# 建立这个文件夹
$ mkdir -p .aspera/connect/etc/
# 复制一份这个文件
$ cp /home/junjun/miniconda3/envs/kingfisher/etc/asperaweb_id_dsa.openssh .aspera/connect/etc/

然后下载就没问题了：

$ kingfisher get -r ERR1739691 -m ena-ascp
10/10/2021 09:02:11 AM INFO: Kingfisher v0.0.1-dev
10/10/2021 09:02:11 AM INFO: Attempting download method ena-ascp ..
10/10/2021 09:02:11 AM INFO: Using aspera ssh key file: $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh
10/10/2021 09:02:11 AM INFO: Querying ENA for FTP paths for ERR1739691..
10/10/2021 09:02:13 AM INFO: Downloading 2 FTP read set(s): ftp.sra.ebi.ac.uk/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_1.fastq.gz, ftp.sra.ebi.ac.uk/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_2.fastq.gz
10/10/2021 09:02:13 AM INFO: Running command: ascp -T -l 300m -P33001 -k 2 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh [e-mail]:/vol1/fastq/ERR173/001/ERR1739691/ERR1739691_1.fastq.gz .

看看后台：

3获取样本信息

接下来获取样本的具体信息，从 SRA Explorer 网址搜索文章的 BioProject 编号 PRJNA634388 ：

然后选择对应样本加入到购物车：

然后下载 full metadata 信息即可：

4下载数据

我们下载的 metadata 信息包括 SRR 号、样本名、fastq 下载地址、aspera 下载地址等等，我们只需要第一列 SRR 号即可：

SRR14765638
SRR14765639
SRR14765640
SRR14765641
SRR14765642
SRR14765643
SRR14765644
SRR14765645
SRR14765646
SRR14765647

写个循环批量下载：

$ mkdir 1.raw-data && cd 1.raw-data

$ batch_download.sh
#!/bin/bash

for i in SRR14765638 SRR14765639 SRR14765640 SRR14765641 SRR14765642 SRR14765643 SRR14765644 SRR14765645 SRR14765646 SRR14765647
do
 kingfisher get -r ${i} -m ena-ascp aws-http prefetch --download-threads 10
 echo "This is ${i} sample download process done !"
done

真是太慢了，还没迅雷快！

干脆用 ascp 下载，我们看看地址，观察其规律看到只有编号在变：

[e-mail]:vol1/fastq/SRR147/038/SRR14765638/SRR14765638_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/038/SRR14765638/SRR14765638_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/039/SRR14765639/SRR14765639_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/039/SRR14765639/SRR14765639_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/040/SRR14765640/SRR14765640_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/040/SRR14765640/SRR14765640_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/041/SRR14765641/SRR14765641_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/041/SRR14765641/SRR14765641_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/042/SRR14765642/SRR14765642_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/042/SRR14765642/SRR14765642_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/043/SRR14765643/SRR14765643_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/043/SRR14765643/SRR14765643_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/044/SRR14765644/SRR14765644_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/044/SRR14765644/SRR14765644_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/045/SRR14765645/SRR14765645_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/045/SRR14765645/SRR14765645_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/046/SRR14765646/SRR14765646_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/046/SRR14765646/SRR14765646_2.fastq.gz
[e-mail]:vol1/fastq/SRR147/047/SRR14765647/SRR14765647_1.fastq.gz
[e-mail]:vol1/fastq/SRR147/047/SRR14765647/SRR14765647_2.fastq.gz

写个循环批量下载：

# 检查
$ for i in {38..47}; \
  do echo "ascp -T -l 300m -P33001 -k 2 \
  -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  [e-mail]:vol1/fastq/SRR147/0${i}/SRR147656${i}/SRR147656${i}_1.fastq.gz . " && \
  echo "This is SRR147656${i} sample R1 download process done !" ;\
  done

# 检查命令没问题就下载
# 下载R1
$ nohup for i in {38..47} ;\
  do ascp -T -l 300m -P33001 -k 2 \
  -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  [e-mail]:vol1/fastq/SRR147/0${i}/SRR147656${i}/SRR147656${i}_1.fastq.gz . && \
  echo "This is SRR147656${i} sample R1 download process done !" ;\
  done &

# 下载R2
$ nohup for i in {38..47} ;\
  do ascp -T -l 300m -P33001 -k 2 \
  -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  [e-mail]:vol1/fastq/SRR147/0${i}/SRR147656${i}/SRR147656${i}_2.fastq.gz . && \
  echo "This is SRR147656${i} sample R2 download process done !" ;\
  done &

更简洁的代码，把地址保存到 sample.txt 里面：

# simple code
$ nohup cat sample.txt |while read id;do ascp -T -l 300m -P33001 -k 2 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh ${id} ./ ; done &

下一节继续...

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,代码已上传至QQ群文件夹，欢迎下载。

群二维码：

老俊俊微信：

知识星球：

所以今天你学习了吗？

欢迎小伙伴留言评论！

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！

 往期回顾 

◀...