wes 7 gost 下载_全外显子测序（WES）分析流程

首先创建一个项目目录wes和存放每个步骤生成文件的子目录raw, clean, align, genome, hg19_VCF

mkdir {raw,clean,align,genome,hg19_VCF}

此流程包含两个raw外显子测序文件：wes.1.fq.gz，wes.2.fq.gz；存放于raw目录

原始数据质控，用fastp

# 项目主目录下运行，在clean目录下生成两个过滤的fq文件
nohup fastp -i raw/wes.1.fq.gz -o clean/wes.1.clean.fq.gz -I raw/wes.2.fq.gz -O clean/wes.2.clean.fq.gz &

下载参考基因组（hg19）文件，存放于genome目录，并建立索引：

for i in $(seq 1 22) X Y M;
do 
nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz & 
done


for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fa;
done


nohup bwa index -a bwtsw -p hg19 hg19.fa  &

比对

# clean目录下运行此命令，在align目录下生成sam文件
nohup bwa mem -t 4 -M -R "@RGtID:lane1tPL:illuminatLB:librarytSM:wes" ../genome/hg19  wes.1.clean.fq.gz  wes.2.clean.fq.gz >../align/wes.sam 2>wes.bwa.align.log &
-t，线程数；
  -M , -M 将 shorter split hits 标记为次优，以兼容 Picard markDuplicates 软件;
  -R 接的是 Read Group的字符串信息，它是用来将比对的read进行分组的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。
  (1) ID，这是Read Group的分组ID，一般设置为测序的lane ID
  (2) PL，指的是所用的测序平台
  (3) SM，样本ID
  (4) LB，测序文库的名字
这些信息设置好之后，在RG字符串中要用制表符（t）将它们分开
samtools view -b -S wes.sam > wes.bam  #sam2bam,便于存储
samtools sort wes.bam -o wes.sorted.bam    #有顺序的排序，便于后面的操作
samtools flagstat wes.sorted.bam > wes.sorted.bam.flagstat    #统计比对信息
1. QC pass的reads的数量为5002344 ，未通过QC的reads数量为0，意味着一共有5002344条reads；
2.标记为secondary的read
3. 嵌合体reads
4.PCR 重复引起的reads
5.  比对到参考基因组上的reads数量；
6.  paired reads数据数量；
7.  read1的数量；
8.  read2 的数量；
9.  正确地匹配到参考序列的reads数量；
10.  一对reads都比对到了参考序列上的数量，但是并不一定比对到同一条染色体上；
11. 一对reads中只有一条与参考序列相匹配的数量； 
12.  一对reads比对到不同染色体的数量；
13.  一对reads比对到不同染色体的且比对质量值大于5的数量。

外显子区域覆盖度

需要生成外显子interval文件，生成这个文件的前提又需要dict文件和外显子bed文件（此处用的是安捷伦外显子bed文件，也可以去UCSC下载）。

gatk CreateSequenceDictionary -R hg19.fa -O hg19.dict
gatk BedToIntervalList -I S31285117_Regions.bed -O Exon.Interval.bed -SD ./genome/hg19.dict
gatk CollectHsMetrics -BI Exon.Interval.bed -TI Exon.Interval.bed -I wes.sorted.bam -O wes.stat.txt