肿瘤外显子数据处理系列教程​（二）质控与去接头

上期的conda安装步骤中没有添加镜像，“conda config --add channels bioconda”，本期的第一节仍旧是数据的下载，对此没有疑问的同学，请直接略过。

因为周一专栏处于6月交接期间，所以大家仍然是可以在头条教程看到新的学徒的数据挖掘教程，然后我这个肿瘤外显子数据处理系列会一直在次页慢慢更新，直到离开为止，非常谢谢大家关注了我一年！

1. 数据下载

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP070662

## 推荐下载NCBI官网最新版本的sra-tools（2.9.6）
cd
mkdir tools
cd tools
curl https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.6-1/sratoolkit.2.9.6-1-ubuntu64.tar.gz --output sratoolkit.2.9.6-1-ubuntu64.tar.gz
tar -zxvf sratoolkit.2.9.6-1-ubuntu64.tar.gz
mkdir bin
cd bin
ls ../sratoolkit.2.9.6-1-ubuntu64/bin/ | grep -e '[A-Za-z]$' | while read id
do
ln -s ../sratoolkit.2.9.6-1-ubuntu64/bin/$id $(basename $id)
done
export PATH='$HOME/tools/bin:$PATH'

## 通过aspera下载sra数据更快
cd
wget 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
export PATH='$HOME/.aspera/connect/bin:$PATH'

source ~/.bashrc

mkdir SRP070662
cd SRP070662
cat SraRunTable.txt | grep "WXS" | cut -f 15 > config
cat SraRunTable.txt | grep "WXS" | cut -f 17 | sed s/_1// >case
paste config case > sra2case.txt

mkdir 0.sra
cd 0.sra
touch download_sra.log

## 官方推荐使用prefetch下载数据
nohup prefetch --option-file ../config  -t fasp -p -c -O . -X 200G >> download_sra.log 2>&1 &

## SRA数据验证
touch validate.out
ls SRR* | grep -v "vdb" | xargs vdb-validate - >> validate.out 2>&1
cat validate.out | grep consistent | wc

4.8G 2月  24  2016 SRR3182418.sra

3.0G 2月  24  2016 SRR3182419.sra

3.1G 2月  24  2016 SRR3182420.sra

3.1G 2月  24  2016 SRR3182421.sra

3.1G 3月   1 10:58 SRR3182422.sra

3.0G 2月  24  2016 SRR3182423.sra

7.8G 3月   1 11:57 SRR3182424.sra

7.8G 3月   4 09:41 SRR3182425.sra

7.7G 2月  24  2016 SRR3182426.sra

8.2G 2月  24  2016 SRR3182427.sra

7.6G 2月  24  2016 SRR3182428.sra

4.1G 2月  24  2016 SRR3182429.sra

3.2G 2月  24  2016 SRR3182430.sra

4.0G 3月   4 09:35 SRR3182431.sra

4.1G 2月  24  2016 SRR3182432.sra

7.6G 2月  24  2016 SRR3182433.sra

6.6G 2月  24  2016 SRR3182434.sra

8.2G 2月  24  2016 SRR3182435.sra

6.0G 2月  24  2016 SRR3182436.sra

6.7G 2月  24  2016 SRR3182437.sra

5.9G 2月  24  2016 SRR3182438.sra

6.2G 3月   1 11:58 SRR3182439.sra

6.0G 2月  24  2016 SRR3182440.sra

7.9G 3月   1 12:03 SRR3182441.sra

6.0G 2月  25  2016 SRR3182442.sra

7.2G 2月  25  2016 SRR3182443.sra

6.8G 2月  25  2016 SRR3182444.sra

5.7G 3月   1 11:42 SRR3182445.sra

7.1G 3月   1 11:42 SRR3182446.sra

2.5G 2月  24  2016 SRR3182447.sra

2. SRA转FQ（fasterq-dump）

## 不建议还使用fastq-dump，速度很慢，但是fasterq-dump的参数不如fastq-dump丰富，没有“-A”参数，所以需要自己修改name，建议最后使用pigz压缩，fasterq-dump和pigz都可以设置多线程

touch sra2fq.log

## sra2fq.sh
cat ../sra2case.txt | while read id
do
   arr=(${id})
   sample=${arr[0]}.sra
   case=${arr[1]}
   echo "time fasterq-dump -p -x -3 -N -P -f --skip-technical -e 16 ${sample} -O . >> sra2fq.log 2>&1"
   echo "sed s/${sample}/${case}/ ${sample}_1.fastq > ${case}_1.fastq "
   echo "pigz -p 16 ${case}_1.fastq"
   echo "sed s/${sample}/${case}/ ${sample}_2.fastq > ${case}_2.fastq "
   echo "pigz -p 16 -f ${case}_2.fastq"
done > sra2fq.sh

nohup bash sra2fq.sh &

5.8G 5月  16 16:57 case1_biorep_A_techrep_1.fastq.gz

5.8G 5月  16 17:04 case1_biorep_A_techrep_2.fastq.gz

4.9G 5月  15 20:46 case1_biorep_B_1.fastq.gz

4.9G 5月  15 20:50 case1_biorep_B_2.fastq.gz

6.1G 5月  15 21:12 case1_biorep_C_1.fastq.gz

6.1G 5月  15 21:18 case1_biorep_C_2.fastq.gz

2.4G 5月  15 18:15 case1_germline_1.fastq.gz

2.4G 5月  15 18:16 case1_germline_2.fastq.gz

5.1G 5月  15 23:13 case1_techrep_2_1.fastq.gz

5.2G 5月  15 23:18 case1_techrep_2_2.fastq.gz

4.5G 5月  15 21:30 case2_biorep_A_1.fastq.gz

4.6G 5月  15 21:34 case2_biorep_A_2.fastq.gz

5.0G 5月  15 21:49 case2_biorep_B_techrep_1.fastq.gz

5.2G 5月  15 21:53 case2_biorep_B_techrep_2.fastq.gz

4.4G 5月  15 22:06 case2_biorep_C_1.fastq.gz

4.5G 5月  15 22:09 case2_biorep_C_2.fastq.gz

2.2G 5月  15 17:45 case2_germline_1.fastq.gz

2.3G 5月  15 17:47 case2_germline_2.fastq.gz

1.4G 5月  18 12:21 case2_germline_2_val_2.fq.gz

4.3G 5月  15 23:32 case2_techrep_2_1.fastq.gz

4.3G 5月  15 23:35 case2_techrep_2_2.fastq.gz

6.0G 5月  15 19:30 case3_biorep_A_1.fastq.gz

6.1G 5月  15 19:36 case3_biorep_A_2.fastq.gz

5.7G 5月  15 19:54 case3_biorep_B_1.fastq.gz

5.6G 5月  15 19:59 case3_biorep_B_2.fastq.gz

3.1G 5月  15 20:07 case3_biorep_C_techrep_1.fastq.gz

3.2G 5月  15 20:09 case3_biorep_C_techrep_2.fastq.gz

2.4G 5月  15 18:00 case3_germline_1.fastq.gz

2.4G 5月  15 18:01 case3_germline_2.fastq.gz

4.4G 5月  16 17:20 case3_techrep_2_1.fastq.gz

4.5G 5月  16 17:23 case3_techrep_2_2.fastq.gz

5.8G 5月  16 16:27 case4_biorep_A_1.fastq.gz

5.8G 5月  16 16:34 case4_biorep_A_2.fastq.gz

5.9G 5月  15 18:36 case4_biorep_B_techrep_1.fastq.gz

5.9G 5月  15 18:42 case4_biorep_B_techrep_2.fastq.gz

5.7G 5月  15 19:03 case4_biorep_C_1.fastq.gz

5.7G 5月  15 19:09 case4_biorep_C_2.fastq.gz

2.4G 5月  15 17:52 case4_germline_1.fastq.gz

2.4G 5月  15 17:54 case4_germline_2.fastq.gz

5.8G 5月  16 17:48 case4_techrep_2_1.fastq.gz

6.0G 5月  16 17:55 case4_techrep_2_2.fastq.gz

4.6G 5月  15 22:24 case5_biorep_A_1.fastq.gz

4.7G 5月  15 22:27 case5_biorep_A_2.fastq.gz

4.6G 5月  16 18:12 case5_biorep_B_techrep_1.fastq.gz

4.7G 5月  16 18:16 case5_biorep_B_techrep_2.fastq.gz

2.0G 5月  16 00:08 case5_biorep_C_1.fastq.gz

1.9G 5月  16 00:09 case5_biorep_C_2.fastq.gz

3.6G 5月  15 17:38 case5_germline_1.fastq.gz

3.7G 5月  15 17:40 case5_germline_2.fastq.gz

5.3G 5月  15 23:56 case5_techrep_2_1.fastq.gz

5.3G 5月  16 00:03 case5_techrep_2_2.fastq.gz

2.4G 5月  16 18:24 case6_biorep_A_techrep_1.fastq.gz

2.5G 5月  16 18:26 case6_biorep_A_techrep_2.fastq.gz

3.1G 5月  15 20:18 case6_biorep_B_1.fastq.gz

3.1G 5月  15 20:20 case6_biorep_B_2.fastq.gz

3.1G 5月  15 20:27 case6_biorep_C_1.fastq.gz

3.1G 5月  15 20:29 case6_biorep_C_2.fastq.gz

2.4G 5月  15 18:07 case6_germline_1.fastq.gz

2.4G 5月  15 18:09 case6_germline_2.fastq.gz

5.2G 5月  15 22:48 case6_techrep_2_1.fastq.gz

5.5G 5月  15 22:53 case6_techrep_2_2.fastq.gz

3. 去接头前质控

## fastqc使用多线程要用--threads

wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.8.zip
unizp fastqc_v0.11.8.zip
ln -s ~/tools/FastQC/fastqc ~/tools/bin/fastqc

mkdir -p 2.qc/{pre, post}
nohup fastqc --outdir ../2.qc/pre --threads 16 *.fastq.gz > ../2.qc/pre/fastqc.log 2>&1 &


## 需要先安装python
## sudo apt install libffi-dev
cd
cd tools/
mkdir python
wget https://www.python.org/ftp/python/3.7.3/Python-3.7.3.tgz
tar zxvf Python-3.7.3.tgz
cd Python-3.7.3
./configure --prefix="/home/llwu/tools/python/"
make
make install
cd ../python/bin
ls * | while read id
do
ln -s ~/tools/python/bin/${id} ~/tools/bin/${id}
done

curl -LOk https://github.com/ewels/MultiQC/archive/master.zip
unzip master.zip
cd MultiQC-master
python3.7 setup.py install

multiqc ./ -n pre_qc_report -p -i " QC REPORT OF SRP070662" -o multiqc


4. 去接头与质控

cd
cd tools
curl https://codeload.github.com/FelixKrueger/TrimGalore/zip/0.6.1 TrimGalore.zip
unzip TrimGalore.zip
ln -s ~/tools/TrimGalore-0.6.1/trim_galore ~/tools/bin/trim_galore

pip3 install --user --upgrade cutadapt
cutadapt --version
## 2.6

cd 1.raw_fq
touch ../3.clean/trimgalore.log

## trim_galore.sh
cat ../tmp | while read id; do
    fq1=${id}_1.fastq.gz
    fq2=${id}_2.fastq.gz
    trim_galore  --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 8 -o ../3.clean  $fq1  $fq2 >> trimgalore.log 2>&1
done

nohup bash trim_galore.sh &

multiqc ./ -n trim_galore_report -p -i " Trim REPORT OF SRP070662" -o multiqc

nohup fastqc --outdir ../2.qc/post --threads 16 *.fq.gz > ../2.qc/post/fastqc.log 2>&1

multiqc ./ -n post_qc_report -p -i " QC REPORT OF SRP070662" -o multiqc


5.0G 5月  18 13:32 case1_biorep_A_techrep_1_val_1.fq.gz

5.1G 5月  18 13:32 case1_biorep_A_techrep_2_val_2.fq.gz

4.4G 5月  17 17:35 case1_biorep_B_1_val_1.fq.gz

4.4G 5月  17 17:35 case1_biorep_B_2_val_2.fq.gz

5.2G 5月  17 17:56 case1_biorep_C_1_val_1.fq.gz

5.3G 5月  17 17:56 case1_biorep_C_2_val_2.fq.gz

2.1G 5月  17 15:23 case1_germline_1_val_1.fq.gz

2.1G 5月  17 15:23 case1_germline_2_val_2.fq.gz

4.3G 5月  17 20:09 case1_techrep_2_1_val_1.fq.gz

4.4G 5月  17 20:09 case1_techrep_2_2_val_2.fq.gz

3.9G 5月  17 18:17 case2_biorep_A_1_val_1.fq.gz

4.0G 5月  17 18:17 case2_biorep_A_2_val_2.fq.gz

4.2G 5月  18 13:11 case2_biorep_B_techrep_1_val_1.fq.gz

4.4G 5月  18 13:11 case2_biorep_B_techrep_2_val_2.fq.gz

3.8G 5月  17 18:34 case2_biorep_C_1_val_1.fq.gz

3.9G 5月  17 18:34 case2_biorep_C_2_val_2.fq.gz

1.3G 5月  17 14:45 case2_germline_1_val_1.fq.gz

1.4G 5月  17 14:45 case2_germline_2_val_2.fq.gz

3.8G 5月  17 20:28 case2_techrep_2_1_val_1.fq.gz

3.8G 5月  17 20:28 case2_techrep_2_2_val_2.fq.gz

5.3G 5月  17 16:26 case3_biorep_A_1_val_1.fq.gz

5.3G 5月  17 16:26 case3_biorep_A_2_val_2.fq.gz

5.0G 5月  17 16:49 case3_biorep_B_1_val_1.fq.gz

5.0G 5月  17 16:49 case3_biorep_B_2_val_2.fq.gz

2.7G 5月  18 12:56 case3_biorep_C_techrep_1_val_1.fq.gz

2.8G 5月  18 12:56 case3_biorep_C_techrep_2_val_2.fq.gz

2.1G 5月  17 15:04 case3_germline_1_val_1.fq.gz

2.1G 5月  17 15:04 case3_germline_2_val_2.fq.gz

3.8G 5月  17 19:31 case3_techrep_2_1_val_1.fq.gz

3.9G 5月  17 19:31 case3_techrep_2_2_val_2.fq.gz

5.1G 5月  17 15:40 case4_biorep_A_1_val_1.fq.gz

5.1G 5月  17 15:40 case4_biorep_A_2_val_2.fq.gz

5.2G 5月  18 12:49 case4_biorep_B_techrep_1_val_1.fq.gz

5.2G 5月  18 12:49 case4_biorep_B_techrep_2_val_2.fq.gz

5.0G 5月  17 16:03 case4_biorep_C_1_val_1.fq.gz

5.0G 5月  17 16:03 case4_biorep_C_2_val_2.fq.gz

2.1G 5月  17 14:54 case4_germline_1_val_1.fq.gz

2.1G 5月  17 14:54 case4_germline_2_val_2.fq.gz

5.0G 5月  17 19:12 case4_techrep_2_1_val_1.fq.gz

5.1G 5月  17 19:12 case4_techrep_2_2_val_2.fq.gz

4.0G 5月  17 18:52 case5_biorep_A_1_val_1.fq.gz

4.0G 5月  17 18:52 case5_biorep_A_2_val_2.fq.gz

3.9G 5月  18 13:20 case5_biorep_B_techrep_1_val_1.fq.gz

3.9G 5月  18 13:20 case5_biorep_B_techrep_2_val_2.fq.gz

1.7G 5月  17 21:01 case5_biorep_C_1_val_1.fq.gz

1.7G 5月  17 21:01 case5_biorep_C_2_val_2.fq.gz

3.1G 5月  17 14:35 case5_germline_1_val_1.fq.gz

3.2G 5月  17 14:35 case5_germline_2_val_2.fq.gz

4.6G 5月  17 20:47 case5_techrep_2_1_val_1.fq.gz

4.6G 5月  17 20:47 case5_techrep_2_2_val_2.fq.gz

2.1G 5月  18 13:01 case6_biorep_A_techrep_1_val_1.fq.gz

2.1G 5月  18 13:01 case6_biorep_A_techrep_2_val_2.fq.gz

2.7G 5月  17 17:06 case6_biorep_B_1_val_1.fq.gz

2.7G 5月  17 17:06 case6_biorep_B_2_val_2.fq.gz

2.7G 5月  17 17:18 case6_biorep_C_1_val_1.fq.gz

2.7G 5月  17 17:18 case6_biorep_C_2_val_2.fq.gz

2.1G 5月  17 15:14 case6_germline_1_val_1.fq.gz

2.1G 5月  17 15:14 case6_germline_2_val_2.fq.gz

4.3G 5月  17 19:50 case6_techrep_2_1_val_1.fq.gz

4.6G 5月  17 19:50 case6_techrep_2_2_val_2.fq.gz

大家可以看到，我这边的分享以干货为主，大概就是人少话不多的代表了吧，不过关于测序数据的质量控制，jimmy老师其实给了我非常信息的指导：

http://www.biotrainee.com/thread-324-1-1.html  （点击下面的阅读原文直达自行学习）

■   ■   ■

• 生信技能树（爆款入门培训课）巡讲第一站-重庆  （已结束）

• 生物信息学全国巡讲之粤港澳大湾区专场 （已结束）

• 生信技能树（爆款入门培训课）巡讲第二站-济南 （已结束）

• 生信技能树（爆款入门培训课）巡讲-千呼万唤进北京（已结束）