3 读文献并下载数据

阅读文献

本文分析的数据来自于文章：Reliability of Whole-Exome Sequencing for Assessing Intratumor Genetic Heterogeneity, Shi et al., 2018, Cell Reports
https://doi.org/10.1016/j.celrep.2018.10.046

原文作者主要是为了对比不同数据处理流程的结果异同与是否可靠，数据处理流程如下：

主要的全外显子组测序有：6 个 case，每个 case 有 5 个样本，分别是一个 Germline，A、B、C 三个 Biological replicate，一个 Technical replicate，一共是 30 个样本。

方法一：下载SRA数据

（不推荐）
NCBI 的 Sequence Read Archive (SRA)，每个项目的 url 格式都是一样的，https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRPXXX 。在文章中，作者给出了 SRA 数据库的 id ：SRP070662

我们可以先去到 sra 数据库：https://www.ncbi.nlm.nih.gov/sra/?term=，检索作者提供的 sra id：SRP070662

点击下图链接

进入该 id 所在的 sra 数据库详细列表
https://www.ncbi.nlm.nih.gov/Traces/study/?query_key=1&WebEnv=NCID_1_75412074_130.14.22.33_5555_1577745764_1047693658_0MetA0_S_HStore&o=acc_s%3Aa

点击 RunInfo Table 和 Accession List 下载两个文件并上传至服务器，存放至 ~/wes_cancer/project 目录下。

cd  ~/wes_cancer/project
ls 
## 0.sra    1.raw_fq    2.clean_fq    3.qc    4.align    5.gatk    6.mutect    7.annotation    8.cnv    9.pyclone    10.signature    SRR_Acc_List.txt    SraRunTable.txt

然后从 SRR_Acc_List.txt 文件中拿到 id 号并用 sra-tools 的 prefetch 命令下载 sra 文件

cat SRR_Acc_List.txt | while read id
do
    prefetch ${id} -O  ./0.sra
done

下载完成后，在 0.sra 目录下，每个 sra 文件以单独的文件夹的形式存放，每一个 sra 文件大小为 2~7G 不等。需要注意的是，还有一些非 sra 格式的文件也会被下载下来，要和 sra 文件放在一起，后面做格式转换的时候程序会自动检索，如果没有放在一起，可能会报错。不过我们这里只关心 sra 文件。

4.8G SRR3182418.sra
3.0G SRR3182419.sra
3.1G SRR3182420.sra
3.1G SRR3182421.sra
......

文件重命名

接下对 sra 文件进行重命名，方法如下：

## 首先提取 sra 文件的 id
cat SraRunTable.txt | grep "WXS" | cut -f 15 > sra
## 然后提取 sra 文件对应的样本 id
cat SraRunTable.txt | grep "WXS" | cut -f 17 | sed s/_1// >config
## 合并，结果如下：两列，第一列为 sra 文件的 id，第二列为重命名后的样本 id
paste sra config > sra2case.txt
cat sra2case.txt
# SRR3182418    case5_germline
# SRR3182419    case2_germline
# SRR3182420    case4_germline
# SRR3182421    case3_germline
# SRR3182422    case6_germline
# SRR3182423    case1_germline
# SRR3182425    case4_biorep_B_techrep
# SRR3182426    case4_biorep_C
# SRR3182427    case3_biorep_A
# SRR3182428    case3_biorep_B
# ......
## 用下面代码重命名
cat sra2case.txt | while read id
do
   arr=(${id})
   sample=${arr[0]}
   case=${arr[1]}
   mv ./0.sra/${sample}/${sample}.sra  ./0.sra/${case}.sra
done
## 最后删掉没用的文件夹
rm -r 0.sra/SRR*
# nohup prefetch --option-file ../sra  -t fasp -p -c -O . -X 50G > download_sra.log 2>&1 &
## vdb-validate
# touch validate.out
# ls SRR* | grep -v "vdb" | xargs vdb-validate - >> validate.out 2>&1
# cat validate.out | grep consistent | wc

通过上面的方法重命名后，得到如下结果：

case1_biorep_A_techrep.sra  
case1_biorep_B.sra            
case1_biorep_C.sra            
case1_germline.sra
case1_techrep.sra
......

格式转换 sra 转 fq

最开始将 sra 转成 fq 用到的命令是 fast-dump，但是速度很慢。新版本的 sra-tools ，可以使用更快的工具做格式转换 fasterq-dump。按照我的代码习惯，需要用到一个 config 配置文件，其实里面就是记录了样本的 id，后面我们会经常用到这个文件。

$ cat config
case1_biorep_A_techrep
case2_biorep_A
case3_biorep_A
case4_biorep_A
case5_biorep_A
case6_biorep_A_techrep
case1_biorep_B
case2_biorep_B_techrep
case3_biorep_B
case4_biorep_B_techrep
case5_biorep_B_techrep
case6_biorep_B
case1_biorep_C
case2_biorep_C
case3_biorep_C_techrep
case4_biorep_C
case5_biorep_C
case6_biorep_C
case1_techrep_2
case2_techrep_2
case3_techrep_2
case4_techrep_2
case5_techrep_2
case6_techrep_2

这里我写了个简单的脚本 sra2fq.sh，然后通过 nohup…& 提交至服务器后台运行。

cd ./0.sra
## fasterq-dump
cat config | while read id
do
    time fasterq-dump -3 -e 16 ${id}.sra -O ../1.raw_fq --outfile ${id}.fastq  >> ../1.raw_fq/${id}_sra2fq.log 2>&1
    time pigz -p 10 -f ../1.raw_fq/${id}_1.fastq >../1.raw_fq/${id}_1.fastq.gz
    time pigz -p 10 -f ../1.raw_fq/${id}_2.fastq >../1.raw_fq/${id}_2.fastq.gz
done 
nohup bash sra2fq.sh &

解释一下上面的脚本，首先通过 fasterq-dump 将 sra 转成 fastq格式，-e 16 指的是调用了 16 个线程，输出到了 1.raw_fq 文件夹中，得到了原始的 fq 文件，因为是双端测序，所以每个样本会有 2 个 fastq 文件${id}_1.fastq 和 ${id}_2.fastq，对应 reads_1 和 reads_2，还有一个垃圾文件 ${id}_fastq。原始的 fastq 文件很大，所以通过 pigz 压缩命令调用了 -p 10 个线程来压缩，可以留意一下压缩后文件的大小。

5.8G  case1_biorep_A_techrep_1.fastq.gz
5.8G  case1_biorep_A_techrep_2.fastq.gz
4.9G  case1_biorep_B_1.fastq.gz
4.9G  case1_biorep_B_2.fastq.gz
6.1G  case1_biorep_C_1.fastq.gz
6.1G  case1_biorep_C_2.fastq.gz
2.4G  case1_germline_1.fastq.gz
2.4G  case1_germline_2.fastq.gz
......

方法二：下载 fastq 数据

上面介绍的 SRA 数据下载比较麻烦，很多读者反映下载遇到问题，评论区更新几种可行的下载方法，这里我就再总结一个新的方法，从 EBI 数据库直接下载 fastq 格式数据，直接越过 SRA，而且不用再做格式转换。

1. 安装 Aspera

在服务器上安装一个 Aspera 软件，可以提高下载速度，方法很简单

cd  ~/wes_cancer/biosoft
wget -c  https://download.asperasoft.com/download/sw/connect/3.8.1/ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
tar -zxvf ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh

安装完成之后，会在家目录生成一个以点号隐藏文件夹 .aspera，为了方便使用，将其添加到环境变量，注意这里是两个大于号>>表示追加，不要打漏了：

echo 'export PATH=~/.aspera/connect/bin/:$PATH' >> ~/.bashrc
source ~/.bashrc

2. 获取下载链接

感谢@城管大队哈队长的提醒，可以到这个链接 https://sra-explorer.info/ 获取我们要下载的数据的 URL，直接在搜索框输入要下载的数据 id 号，然后选择要下载的数据即可，这里我仅仅下载外显子测序数据：

然后给出几种下载方式，这里我们直接选择 fastq 格式，用 Aspera 下载，网页就会自动生成下载的方式，直接copy或者download到服务器上就好，这里它还很个性化的帮我们重命名了，但是我后面还会再改名的：

如果你没法通过 SRA-explorer 网页获取下载地址，就直接去 ENA 数据库搜索 id 号然后找到下载地址，就可以通过 ascp 命令来下载了。

3. 下载数据

这里我是从上面网页 copy 下载方式的，然后在服务器上粘贴一下，写入一个脚本中：

$ cd ~/wes_cancer/project/1.raw_fq
$ cat >fq_download.sh
#!/usr/bin/env bash
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/008/SRR3182418/SRR3182418.fastq.gz . && mv SRR3182418.fastq.gz SRR3182418_exome_sequencing_of_case5_germline.fastq.gz
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/008/SRR3182418/SRR3182418_1.fastq.gz . && mv SRR3182418_1.fastq.gz SRR3182418_exome_sequencing_of_case5_germline_1.fastq.gz
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/008/SRR3182418/SRR3182418_2.fastq.gz . && mv SRR3182418_2.fastq.gz SRR3182418_exome_sequencing_of_case5_germline_2.fastq.gz
......
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/006/SRR3182446/SRR3182446_2.fastq.gz . && mv SRR3182446_2.fastq.gz SRR3182446_exome_sequencing_of_case5_techrep_2_2.fastq.gz
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/007/SRR3182447/SRR3182447.fastq.gz . && mv SRR3182447.fastq.gz SRR3182447_exome_sequencing_of_case5_biorep_C.fastq.gz
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/007/SRR3182447/SRR3182447_1.fastq.gz . && mv SRR3182447_1.fastq.gz SRR3182447_exome_sequencing_of_case5_biorep_C_1.fastq.gz
ascp -QT -l 300m -P33001 -i $HOME/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR318/007/SRR3182447/SRR3182447_2.fastq.gz . && mv SRR3182447_2.fastq.gz SRR3182447_exome_sequencing_of_case5_biorep_C_2.fastq.gz
^C

然后直接运行这个脚本就好：

bash fq_download.sh

至于速度，大家自己看吧：

当然，这么大的数据，这样下载占用了命令行，所以，通常我们都是挂在后台下载的：

nohup  bash  fq_download.sh  &

下载完成后，每个样本会有 3 个 fastq 文件，分别对应 *_1.fastq.gz 和 *_2.fastq.gz ，以及一个只有几百K没什么用的 fastq 文件，比如下面这个样本：

442K SRR3182418_exome_sequencing_of_case5_germline.fastq.gz
3.3G SRR3182418_exome_sequencing_of_case5_germline_1.fastq.gz
3.4G SRR3182418_exome_sequencing_of_case5_germline_2.fastq.gz

通常这个多余的 fastq 文件我都不关心，直接删掉，一个个删掉也有点麻烦，当然如果你的 Linux 基础比较好的话，是可以用命令来删除的：

find ./*gz -size 1M | xargs rm

4. 重命名

虽然上面网页给的下载方式会顺便帮我们重命名，但是还是不符合我的命名习惯，因此，这里我再重新命名一次，通过观察文件名，可以发现每个文件的前 31 个字符都是可以去掉的，我们取后半部分就好，后半部分长短不一，我就直接取 100 个字符长度吧，管它呢。一句简单的代码即可实现：

ls *gz | while read id; do mv ${id} ${id:31:100}; done

结果是：

case1_biorep_A_techrep_1_1.fastq.gz  
case1_biorep_A_techrep_1_2.fastq.gz  
case1_biorep_B_1.fastq.gz            
case1_biorep_B_2.fastq.gz        
......

注意一下，这里的技术重复样本，和我后面文章的命名方式稍有不同，如果读者是用这种方法下载数据，请在将后文经常用到的 config 文件进行一定的修改，至于怎么修改，就自己摸索一下吧：

case1_biorep_A_techrep_1_1.fastq.gz  
case1_biorep_A_techrep_1_2.fastq.gz  
case1_techrep_2_1.fastq.gz           
case1_techrep_2_2.fastq.gz