[TOC]
0. 背景和准备
WES测序原理和过程
实验流程
实验流程
数据分析流程
数据分析流程
标准信息分析
- 去污染、去低质量等过滤处理
- 数据产出的统计
- 外显子区域测序深度的直方分布图
- SNPs的检测
- SNPs的注释
- 插入和缺失（Indels）的检测
- 插入和缺失（Indels）的注释
个性化信息分析
- 氨基酸替代的预测分析
- 群体SNP的获取和等位基因频率评估
- Mendelian disorder analysis 孟德尔遗传疾病分析
- 基于新一代测序技术的全基因组关联（NGS-GWAS）分析
- 外显子融合分析
- 基因融合
文献和综述阅读
Exome sequencing covers >98% of mutations identified on targeted next generation sequencing panels
微信推文：(腾讯不时更换推文网址，不能保证链接有效)
基因检测与解读:全外显子组测序在临床和科研中的现状、需求和选择攻略
基因谷 :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格式
参考阅读：
每4行一条read，一个read包含内容：
id
sequence
+
quality score
quality score的转换：
Phred + 33对应ascii值，转换成ascii码正好对应一个字符，只占一个位置，可与碱基一一对应。
还有一种是phred + 64的换算，但采用这一体系的都是两年前老数据，如有需要可读相关资料。一种简单的判断方法是，如果每个reads的第4行的质量评分字符绝大部分是大写字母，则应是phred + 33。
1541065388751.png
表示该碱基测序准确的概率
下表表示换算的例子：
| 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
网友评论