WES

2018-11-04 本文已影响104人 Stone_Stan4d

[TOC]

0. 背景和准备

WES测序原理和过程

全外显子组测序-一*简介

实验流程

数据分析流程

标准信息分析

去污染、去低质量等过滤处理
数据产出的统计
外显子区域测序深度的直方分布图
SNPs的检测
SNPs的注释
插入和缺失（Indels）的检测
插入和缺失（Indels）的注释

个性化信息分析

氨基酸替代的预测分析
群体SNP的获取和等位基因频率评估
Mendelian disorder analysis 孟德尔遗传疾病分析
基于新一代测序技术的全基因组关联（NGS-GWAS）分析
外显子融合分析
基因融合

文献和综述阅读

Comparison of Exome and Genome Sequencing Technologies for the Complete Capture of Protein‐Coding Regions

Exome sequencing covers >98% of mutations identified on targeted next generation sequencing panels

微信推文：(腾讯不时更换推文网址，不能保证链接有效)

生信技能树:一个WES实战-生信菜鸟团博客2周年精选文章集

基因检测与解读:全外显子组测序在临床和科研中的现状、需求和选择攻略

基因谷 :1168名患者证实！全外显子组测序可大幅改善患者临床治疗选择

目标和流程

WES（一）测序质量控制
WES（二）snp-calling
WES（三）snp-filter
WES（四）不同个体的比较
WES（五）不同软件比较
WES（六）用annovar注释
WES（七）看de novo变异情况

软件和数据准备

软件安装

登陆服务器后，在~目录下创建biosoft文件夹，进入该文件夹，进行后续操作。

miniconda安装

wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
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

查看当前conda环境

conda info --envs

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

source activate wes

安装 sra-tools软件

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

后续有软件需要时安装。

cd ~
mkdir /vip/cyshi/WES
ln -s /vip/cyshi/WES WES 
#因为~目录下空间不够，所以在/vip下创建cyshhi/WES的目录，并建软链接到家目录下
#如果~目录硬盘足够，可以忽略这一步，直接到下一步。

创建一个文件夹WES：

#mkdir ~/WES  #如果上一个代码框建链接的代码未执行，则要运行此行代码；反之，略过。
cd ~/WES
mkdir {raw, clean, qc, align,mutation}
cd ~/WES/raw  #进入~/WES/raw目录

数据下载

这里，先生成一个srrList.txt

SRR1139999  NPC10F-T
SRR1140007  NPC10F-N
SRR1139956  NPC15F-T
SRR1139958  NPC15F-N
SRR1139966  NPC29F-T
SRR1139973  NPC29F-N
SRR1140015  NPC34F-T
SRR1140023  NPC34F-N
SRR1140044  NPC37F-T
SRR1140045  NPC37F-N

然后写个小脚本, 取名为1.down2fq.sh：

#! /bin/bash
#$1 传入每行第一列是要获取的srr号， 第二列是样本名的文本文件名
#$2 传入*fastq.gz输出的文件夹
cat $1 | cut -f 1 | while read srr
do
        nohup prefetch $srr &

done

cat $1 |while read line
do
        array=($line)
        name=${array[0]}
        sample=${array[1]}
        nohup fastq-dump --gzip --split-3 -A $sample ~/ncbi/public/sra/${name}.sra -O $2 &
done

在~/WES文件夹下运行脚本：

nohup bash 1.down2fq.sh srrList.txt ./raw &

下载的数据是双端测序，共5个样本，每个样本分为诊断为鼻咽癌（NPC）的组织和正常组织。输出结果如下：

1541174485735.png

)

1. 质量控制

fastq格式

参考阅读：

fastq格式详解

phred质量评分

每4行一条read，一个read包含内容：

id
sequence
+
quality score

quality score的转换：

$phred= -10lg(1 - p)$

$Phred33 = phred + 33$

Phred + 33对应ascii值，转换成ascii码正好对应一个字符，只占一个位置，可与碱基一一对应。

还有一种是phred + 64的换算，但采用这一体系的都是两年前老数据，如有需要可读相关资料。一种简单的判断方法是，如果每个reads的第4行的质量评分字符绝大部分是大写字母，则应是phred + 33。

$ascii(A) = 65$

$ascii(a) = 97$

1541065388751.png

$p$ 表示该碱基测序准确的概率

ASCII码对照表

下表表示换算的例子：

p值	1-p值	phred质量评分	ascii值	ascii码	含义(错误概率)
99%	0.01	20	53	5	1 error in 100 base
99.9%	0.001	30	63	？	1 error in 1000 base
99.99%	0.0001	40	73	I	1 error in 1000 base
99.999%	0.00001	50	83	S	1 error in 10000 base

Quality Control

质控目的：

了解数据质量、大小、
过滤数据
- 去除接头(画图示)
- 去除低质量碱基(-q 25)
- 最大允许错误率(默认-e 0.1)
- 去除<36的reads(--length 36)
- 切除双端的overlap>3的碱基
- 去除reads以对为单位(--paired)

软件：

fastqc
cutadapt
Trim Galore

###1.data statistics 
qc="~/WES/qc"
raw="~/WES/raw"
nohup find $raw -name *.gz |xargs fastqc -t 10 -o $qc/ &
multiqc *.zip
###2.filter data
trim_galore --phred33  -q 25 -e 0.1  --length 36 --stringency 3 --paired -o ./ $qc  ${raw}/NPC10F-T_1.fastq.gz ${raw}/NPC10F-T_2.fastq.gz

multiqc生成的文件如下：

1541209728318.png

打开网页：

1541209771201.png

批量质控

脚本名字，2.qc2val.sh，内容：

#! /bin/bash
cat $1 |while read id
do 
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
echo $fq1 $fq2
trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ./qc  $qc  $fq1 $fq2
done

这里要传递一个文件列表的参数，每一行为一个样本的配对的双端测序fq文件，此时，在~/WES路径下，可如此生成：

ls ./raw/*fastq.gz | sort | paste - - | > fqPairedList.txt

运行脚本：

bash 2.qc2val.sh  fqPairedList.txt

异常的来源

2 比对

比对目的：

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

比对策略：

index先建立索引
mem算法(maximal exact matches)
跨外显子比对

软件：

bwa(Burrows-Wheeler Aligner )
hisat2

参考基因组：

hg38.fa

人类基因组fa文件下载地址：

http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/

1541071773964.png

相关参考：个人基因组比对及其变异分析

### 1 make.index 可用别人建好的
#### bwa index
cd  ~/reference/index/bwa/hg38
ln -s /public/biosoft/GATK/resources/bundle/hg38/Homo_sa piens_assembly38.fasta ./hg38.fa 
bwa index -a ./hg38.fa
#### hisat2 index
#下载网址https://ccb.jhu.edu/software/hisat2/index.shtml
#### subjunc index
#### bowtie2 index

回帖到基因组

# 样本id 
id=NPC10F-N
# 不同比对策略 

# hisat2-build建索引 
hisat2 -p 10 -x ~/reference/index/hisat/hg38/genome -1 ./qc/${id}_1_val_1.fq.gz -2 ./qc/${id}_2_val_2.fq.gz -S ${id}.hisat.sam 2>./align/${id}.sam.log 

# subread-buildindex 建索引 
subjunc -T 5 -i ~reference/index/subread/hg38 -r ./qc/${id}_1_val_1.fq.gz -R ./qc/${id}_2_val_2.fq.gz -o ./align/${i d}.subjunc.sam
#bowtie2-build建索引 
bowtie2 -p 10 -x ~/reference/index/bowtie/hg38 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ./align/${id}.bowtie.sam 
# bwa index建索引 
bwa mem -t 5 -M ~/reference/index/bwa/hg38 ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz >./align/${id}.bw a.sam

bam排序

cat $1 | while read id
do
        array=($id)   #括号很重要
        fq1=${array[0]}
        fq2=${array[1]}
        sample=${array[2]}
        echo $fq1 $fq2 $sample
        hisat2 -p 10 -x ~/reference/index/hisat2/hg3
8/genome  -1 ./raw/$fq1 -2 ./raw/fq2
-S ./align/${sample}.hisat.sam  -
        echo "$sample map and sort done"
done

WES