2020-07-09 WGS数据分析处理流程

1.拿到数据后先检查数据是否完整。用md5sum命令。

#生成md5文件
ls KPGP*| while read KPGP; do echo $KPGP;md5sum ${KPGP} >> ${KPGP}.md5; done

#检查完整性，全部显示OK即可
md5sum -c *.md5

2.对数据进行质检。

#质检
nohup fastqc -o /data/XXXX/WGS/01fastqc -t 10 *.fq.gz &

#multiqc合并质检报告查看
nohup multiqc * -o /data/XXXX/WGS/01fastqc/multiqc &
#若不合格还要进行用cutadaptor去接头等操作。此处合格，则不赘述。

3.与参考基因组进行比对。
为什么要比对？我们已经知道NGS测序下来的短序列（read）存储于FASTQ文件里面。虽然它们原本都来自于有序的基因组，但在经过DNA建库和测序之后，文件中不同read之间的前后顺序关系就已经全部丢失了。因此，FASTQ文件中紧挨着的两条read之间没有任何位置关系，它们都是随机来自于原本基因组中某个位置的短序列而已。
因此，我们需要先把这一大堆的短序列捋顺，一个个去跟该物种的参考基因组比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。这也是核心流程真正意义上的第一步，只有完成了这个序列比对我们才有下一步的数据分析。
我们这里将用于流程构建的BWA就是其中最优秀的一个，它将BW(Burrows-Wheeler)压缩算法和后缀树相结合，能够让我们以较小的时间和空间代价，获得准确的序列比对结果。

#bwa建立索引
bwa index -a bwtsw chrom.37.fa 

#bwa men
vi bwa.sh

for(( i=1 ; i<=6 ; i++ ))
do
  bwa mem -t 4 /data/XXXX/WGS/02hg19/chrom.37.fa \
   /data/XXXX/WGS/KPGP-00001_L${i}_R1.fq.gz \
   /data/XXXX/WGS/KPGP-00001_L${i}_R2.fq.gz > /data/XXXXX/WGS/03BWA/L${i}.sam 
done

nohup sh bwa.sh &

4.sam转为bam，并sort好。

#samtools sam2bam+sort
for(( i=1 ; i<=6 ; i++ ))
do
  samtools view L${i}.sam>L${i}.bam  
done

for(( i=1 ; i<=6 ; i++ ))
do
  samtools sort L${i}.sam L${i}.sort
done

为什么需要排序？FASTQ文件里面这些被测序下来的read是随机分布于基因组上面的，第一步的比对是按照FASTQ文件的顺序把read逐一定位到参考基因组上之后，随即就输出了，它不会也不可能在这一步里面能够自动识别比对位置的先后位置重排比对结果。因此，比对后得到的结果文件中，每一条记录之间位置的先后顺序是乱的，我们后续去重复等步骤都需要在比对记录按照顺序从小到大排序下来才能进行。
5.标记PCR重复。
为什么要去除这个duplication？主要是因为在call snp的时候，如果某个变异位点的变异碱基都是来自于PCR重复，而我们却认为它深度足够判断是真的变异位点，这个结论其实有很大可能是假阳性。

#安装好gatk
for(( i=1; i<7; i++ ))
do
  gatk MarkDuplicates -I L${i}.sort.bam -O /data/lanyunzhou/WGS/04gatk.mark/L${i}.sort.markdup.bam -M /data/lanyunzhou/WGS/04gatk.mark/L${i}.markdup_metrics.txt
done

#在bwa的时候没有添加，因此在这一步添加read group
for(( i=1; i<7; i++ ))
do
  java -jar /data/XXXXX/software/software/picard.jar AddOrReplaceReadGroups I=L${i}.sort.markdup.bam  O=L${i}.sort.markdup.addrg.bam  ID=group1 LB= lib1 PL=illumina PU=unit1 SM=sample1
done

#建立新的bam文件的索引
for(( i=1; i<7; i++ ))
do
  samtools index L${i}.sort.markdup.addrg.bam
done

Read Group，是一个非常重要的信息，以@RG开头，它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的。在Read Group中，有如下几个信息非常重要：
(1)ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
(2) PL，指的是所用的测序平台，这个信息不要随便写！特别是当我们需要使用GATK进行后续分析的时候，更是如此！这是一个很多新手都容易忽视的一个地方，在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台，
(3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本。
(4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；
除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后，在RG字符串中要用制表符（\t）将它们分开。

bwa mem -t 4 -R ‘@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name’ /path/to/human.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam

6.初步variants calling。（GATK best software for snp+indel）

#生成参考序列索引
samtools faidx /data/XXXX/WGS/02hg19/chrom.37.fa 

##建参考序列dict
java -jar /data/XXXXX/software/software/picard.jar CreateSequenceDictionary R=/data/XXXXX/WGS/02hg19/chrom.37.fa O=/data/XXXXX/WGS/02hg19//chrom.37.dict

#用gatk的时候才可省去局部重比对。用其他软件不行
#初步calling
for(( i=1; i<7; i++ ))
do
  gatk HaplotypeCaller \
  -R /data/XXXXX/02hg19/chrom.37.fa \
  -I /data/XXXXX/WGS/04gatk.mark/L${i}.sort.markdup.addrg.bam \
  -O /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf
done

7.过滤质量低的变异。

for(( i=1; i<6; i++ ))
do
  bgzip -f L${i}.vcf
done

for(( i=1; i<6; i++ ))
do
  tabix -p vcf L${i}.vcf.gz
done

#SNP
for(( i=1; i<7; i++ )) 
do
  gatk VariantRecalibrator \
  -R /data/XXXXX/WGS/02hg19/chrom.37.fa \
  -V /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf.gz \
  -resource:hapmap,known=false,training=true,truth=true,prior=15.0 /data/XXXXX/WGS/06VQSR.file/hapmap_3.3.hg19.sites.vcf \
  -resource:omini,known=false,training=true,truth=false,prior=12.0 /data/XXXXX/WGS/06VQSR.file/1000G_omni2.5.hg19.sites.vcf \
  -resource:1000G,known=false,training=true,truth=false,prior=10.0 /data/XXXXX/WGS/06VQSR.file/1000G_phase1.snps.high_confidence.hg19.sites.vcf \
  -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /data/XXXXX/WGS/06VQSR.file/dbsnp_138.hg19.vcf \
  -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
  -mode SNP \
  -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
  --rscript-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.plots.R \
  --tranches-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.tranches \
  -O /data/XXXXX/WGS/07gatk.snp/L${i}.snps.recal
done

for(( i=1; i<7; i++ ))
do
  gatk ApplyVQSR \
  -R /data/XXXXX/WGS/02hg19/chrom.37.fa \
  -V /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf.gz \
  -ts-filter-level 99.0 \
  --tranches-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.tranches \
  -recal-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.recal \
  -mode SNP \
  -O /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz
done

#INDEL
for(( i=1; i<7; i++ ))
do
  gatk VariantRecalibrator \
  -R /data/XXXXX/WGS/02hg19/chrom.37.fa \
  -O /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.recal \
  -V /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz \
  -resource:mills,known=true,training=true,truth=true,prior=12.0 /data/XXXXX/WGS/06VQSR.file/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
  -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /data/XXXXX/WGS/06VQSR.file/dbsnp_138.hg19.vcf \
  -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
  -mode INDEL \
  --max-gaussians 6 \
  --rscript-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.plots.R \
  --tranches-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.tranches
done

for(( i=1; i<7; i++ ))
do
  gatk ApplyVQSR \
  -R /data/XXXXX/WGS/02hg19/chrom.37.fa \
  -V /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz \
  -ts-filter-level 99.0 \
  --tranches-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.tranches \
  -recal-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.recal \
  -mode INDEL \
  -O /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indel.VQSR.vcf.gz
done

#解压vcf.gz
bgzip -d *.vcf.gz

8.snpEff对变异进行注释。

#snpeff (到snpeff.jar那个文件夹)
java -jar $snpeff GRCh37.75 /data/XXXXX/WGS/07gatk.snp/L1.snps.indel.VQSR.vcf > /data/XXXXX/WGS/09snpeff.note/L1.snp.snpeff.vcf

9.突变频谱分析（6种）

cat L1.snpeff.vcf |grep -v INDEL |grep -v "^#" |cut -f 1-5 >L1.snps.txt
cat L1.snps.txt |grep -v INDEL |grep -v "^#" |cut -f 4,5|sort |uniq -c |grep -v ",">L1.snps2.txt

dat <- data.frame(type=c('C>A(G>T)','C>T(G>A)','C>G(G>C)','T>A(A>T)','T>G(A>C)','T>C(A>G)'),
                  counts=c(205713+208287,495660+496650,130847+130034,113079+112752,132902+131814,497236+497094))
library(ggplot2)
p=ggplot(dat,aes(x=type,y=counts,fill=type))+geom_bar(stat="identity")
print(p)

10.祖源分析
祖源分析是用千人基因组数据挑位点和自己的基因组数据位点进行PCA分析。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4784403/
人类的Y染色体拥有约0.58亿个碱基对（DNA基本结构），约占人类男性体细胞中DNA的2%。
人类Y染色体上有86个基因，这些基因只编码了23种不同的蛋白质。只有拥有Y染色体才能可能继承的性状被称为雄性性状。
人类的Y染色体除了在端粒上的拟常染色体区的少部分片段（只占有染色体长度约5%）能与相应的X染色体发生重组，
其外都不能发生重组。
这些片区是由原本X染色体与Y染色体同源的片段遗留下来的。
Y染色体中不能发生重组的其他部分被称为“NRY区”（non-recombining region，非重组区）。
这个区域中的单核苷酸多态性被用于父系祖先的追溯。
11.SV(CNS)

#delly
delly call -x hg19.excl -o delly.bcf -g hg19.fa input.bam
bcftools view delly.bcf > delly.vcf

对检测到的SV进行genomic feature的注释

Bioconductor 的intansv包

参考资料：WGS数据分析处理流程