实践WES

WES流程图.PNG

测序数据的输入文件为raw fastq，质控软件为fastqc 和trim_galore，通过参数设置，得到结果文件clean fastq。

比对。输入文件为clean/HQ fastq，比对软件使用BWA（不止这一种），比对算法MEM，结果文件为raw sam，然后转成bam文件。

找变异。IGV可视化bam，GATK过滤bam，GATK找变异得到vcf

注释。输入文件为HQ vcf，比对软件Annovar，使用参考数据库，得到结果文件即与功能相关的变异信息。

后续进行个性化分析。

实战阶段

软件安装

wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-4.5.11-Linux-x86_64.sh
bash Miniconda2-4.5.11-Linux-x86_64.sh

更改环境变量

echo 'export PATH=/home/jianmingzeng/biosoft/myBin/bin:$PATH' >>~/.bashrc 
这个/home/jianmingzeng/biosoft/myBin/bin是个例子，得要找到自己安装了的那个软件的位置。怎么找到？以安装GTK软件为例，先 ./gatk --help 能出来，然后 pwd 出来路径，例如 /home/vip18/GATK4/gatk-4.0.11.0  然后将这个添加到环境变量里。最后这么使之生效：
source ~/.bashrc

这事也可以这么干：
进入vim编辑，
vim ~/.bashrc
写进去，
export PATH="/home/jianmingzeng/biosoft/myBin/bin:$PATH"
使之生效，
source ~/.bashrc

配置镜像

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

创建名为wes的软件安装环境

conda create -n wes python=2 #这里同时安装了python 2版本的软件。也可以后续安装它。

查看当前conda环境

conda info --envs

激活/进入conda的wes环境（避免每次用-n wes）

source activate wes

安装sra-tools软件

conda search sra-tools   # fastq-dump --help
conda install -y sra-tools   # done正确安装，且能调出软件help

继续安装软件

conda install -y fastqc trim-galore star bwa hisat2 bowtie2 subread tophat htseq bedtools deeptools salmon cutadapt multiqc # 注意核查每个软件是否能够被调用

数据下载

找到SRA数据

https://www.ncbi.nlm.nih.gov/sra?term=SRP077971 [可修改SRP号]

找到SRA数据集，找到_List.txt

https://www.ncbi.nlm.nih.gov/Traces/study/?WebEnv=NCID_1_5569658_130.14.22.76_5555_1540039
270_2975166503_0MetA0_S_HStore&query_key=1

获得fastq

source activate wes # step1 激活环境

cat SRR_Acc_List.txt | while read id; do (prefetch ${id});done  # step2获得SRR_Acc_List.txt 或 观察数据编号的规律下载

ls ./1.SRA/*sra | while read id; do ( nohup fastq-dump --gzip --split-3 -O ./ ${id} ); done # step3 SRA数据转成fastq

fastq 格式

fastq格式.PNG

现在用预先下载好的文件做练习，例如
7E5239.L1_1.fastq.gz
7E5241.L1_2.fastq.gz

ll ./*/*fastq.gz  # 可以搜索找到gz结尾的文件。

fastqc软件对练习数据进行质控

fastqc -t 3 ./7E5239.L1_1.fastq.gz  ./7E5241.L1_2.fastq.gz

质控结果会反映在每一个.html文件里，样品较多看起来不方便，可使用multiqc整合在一个.html文件里。multiqc不止能整合fastqc软件运行结果，也可以整合其他软件运行结果。更多信息需官网查询。

运行fastqc时候可能会报错，例如出现Permission denied 字样，表明没有相应权限，这说明该软件结果的输出不在自己的目录中。查询帮助即 fastqc --help 查找出入参数用法，将输出路径更改到自己的目录中来，即定义输出路径。例如：
fastqc /home/qmcui/test_boy/1.fastq_qc/*.fastq.gz -o ./project/wes/

所有样品的fastqc运行成功后再运行multiqc。可先multiqc --help，查看帮助文档。如下运行：
multiqc ./

查看结果，得知reads大小，作为参考，以备运行trim_galore过滤数据，实际该步骤通常不需要自己做，过程较为复杂，通常是测序公司做好即可，以下trim_galore的运行仅为了解。

trim_galore --phred33 -q 25 -e 0.1--length 36 --stringency 3--paired -o ./ 
./7E5241.L1_1.fastq.gz ./7E5241.L1_2.fastq.gz

批量质控。多个样本时候需要批量质控，写脚本文件后缀为.sh，例如取名为trim_galore.sh

cat > 自己的目录 ./trim_galore.sh   # 例如，$ cat /home/qmcui/test_boy/1.fastq_qc/trim_galore.sh 

# 以下为写入文本的内容
cat config|while read id
do 
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
    echo "################################################" # 设定这个是为调试程序。另外，这一长串比较显眼，容易在屏幕上找到。这并不重要。
    trim_galore -q 25 --phred33 --length 36 \
    --stringency 3 --paired -o ./  $fq1 $fq2
    echo "OK"
done

# crtl +c结束

# 在运行之前，先cat 一下刚才写的文件，再查看一下所写内容正确与否，确认无误后再运行。

# 提交后台运行
nohup bash trim_galore.sh & v# 会返回一个后台任务id号

jobs -l   #可以看到刚才返回的后台任务id已经在运行

Mapping

比对目的：
• 将打断测序的reads比回参考基因组
• 得到比对结果sam文本，用于后续分析

比对策略：
• index先建立索引
• mem算法(maximal exact matches)

软件：
• bwa(Burrows-Wheeler Aligner )

参考基因组：
• hg38.fa

Sort

sort目的：
• 用于后续软件分析
• 提取bam序列

比对策略：

一步到位生成sort.bam （bam文件排序）
• bwa生成sam转bam,然后再排序。
• bwa比对时直接利用管道符‘|’直接生成排序bam

软件：
• samtools(sort参数)

一步到位生成sort.bam

bwa mem -t 5 -R "@RG\tlD:$sample\tSM:$sample\tLB:WGS\tPL:illumina" $INDEX $fq1 $fa2 |samtools sort -@ 5 -o $sample.sort.bam -
 

bam统计

ls *.bam|while read id;do(samtools flagstat $id > $(basename$id".bam").stat);done

Mark PCR重复(Mark Duplicates)

Mark Duplicates目的：
• 标记/删除PCR重复的reads
• 为后续call变异位点增加可信度,去掉假阳性

使用软件GATK4

建立路径

mkdir 4.gatk_markdup
cd 4.gatk_markdup

mkdup bam donot delete

sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \
-I ../3.sort_bam/$sample.sort.bam \
-O ${sample}_marked.bam \
-M ${sample}.metrics \
1>${sample}_log.mark 2>&1

找变异——variant call之gatk多样本

参数：
• sample赋值
• -I 输入文件
• -O 输出文件
• 1> 重定向标准输出
• 2> 重定向错误输出到1的文件内(目的：记录程序运行日志)

目的：
• 得到vcf前矫正的bam步骤

软件：
• GATK 子命令FixMateInformation
• samtools 子命令index 建索引

step 1

### 建立路径 ###
mkdir 6.gatk_bam
cd 6.gatk_bam

### 找变异 ###
sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation \
-I ${sample}_marked.bam \
-O ${sample}_marked_fixed.bam \
-SO coordinate \
1>${sample}_log.fix 2>&1
samtools index ${sample}_marked_fixed.bam

step 2

目的：
• 得到vcf前矫正的bam步骤

参数：
• sample赋值样本ID
• ref赋值参考基因组
• snp赋值dbsnp数据库(同indel)
• -I 输入文本
• -O 输出文本(承前接后)

软件：
• GATK 子命令BaseRecalibrator

### 进入路径 ###
cd 6.gatk_bam

### 找变异，先赋值 ###

sample="7E5241"
ref="/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly
38.fasta"
snp="/public/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz"
indel="/public/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold
_standard.indels.hg38.vcf.gz"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator \
-R $ref \
-I ${sample}_marked_fixed.bam \
--known-sites $snp \
--known-sites $indel \
-O ${sample}_recal.table

step 3

参数：
• sample赋值样本ID
• ref赋值参考基因组
• snp赋值dbsnp数据库
• -I 输入文本
• -O 输出文本(承前接后)
目的：
• 得到vcf前矫正的bam步骤

软件——GATK 子命令
• ApplyBQSR
• HaplotypeCaller

### 矫正bam，先赋值 ####
sample="7E5241"
ref="/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly
38.fasta"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR \
-R $ref \
-I ${sample}_marked_fixed.bam \
-bqsr ${sample}_recal.table \
-O ${sample}_bqsr.bam \
1>${sample}_log.ApplyBQSR 2>&1

step 4

### calling variation得到raw vcf ###
mkdir gatk_single_vcf && cd gatk_single_vcf
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
-R $ref \
-I ${sample}_bqsr.bam \
--dbsnp $snp \
-O ${sample}_raw.vcf \
1>${sample}_log.HC 2>&1

raw vcf filter

步骤较多，暂时跳过

注释 （使用annovar）

# 建立路径（没有的话，要先建一个）
# mkdir annovar

~qmcui/software/annovar/annovar/convert2annovar.pl \
-format vcf4old \      # -filter  # -allsample  #另有许多参数，详见官网
${sample}_filtered.vcf > ${sample}.annovar

# 运行完成之后得到.anno.variant_function的文件，可以对这个文件进行进一步的操作，例如：

cat  7E5241.anno.variant_function

less -S  7E5241.anno.variant_function|cut -f 2 | sort|uniq -c #  统计第二列

less -S  7E5241.anno.variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|cut -f 2 | sort|uniq -c

less -S  7E5241.anno.variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|less -S >##.txt

less -S 7E5240_last_annovar.exonic_variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|awk '$4=="chr1" && $5 >0 && $6 <208559422 {print $0}'|less -S
`