利用重测序数据Call SNP,并尝试构建遗传图谱

#构建ref genome index file 

bwa index genome.fasta genome

---

vim sam_file.sh

#!/usr/bin/env bash

for i in {s1,d1,fuben,muben}

do

#将重测序的clean data 回帖到参考基因组上

bwa mem -t 28 -M ./data/genome.fasta ./data/${i}_1.fq ./data/${i}_2.fq >./${i}/${i}.sam

#转换为bam文件

samtools view -@ 28 -S -b  ./${i}/${i}.sam >./${i}/${i}.bam

#统计比对结果

samtools flagstat ./${i}/${i}.bam >./${i}/${i}.log

#进行sort

samtools sort -@ 28 ./${i}/${i}.bam ./${i}/${i}.sort

done

---

nohup bash sam_file.sh &

#使用GATK4 call variants

vim GATK.sh

#!/usr/bin/bash

for i in {s1,d1,fuben,muben}

do

#将sort后的bam文件加上head信息

java -jar /biosoftware/picard-tools-1.134/picard.jar AddOrReplaceReadGroups I=${i}/${i}.sort.bam O=${i}/${i}_picard.bam SO=coordinate ID=${i} LB=${i} SM=${i} PL=illumina PU=run

#建立 bam索引和参考基因组的库文件

samtools index ${i}/${i}_picard.bam

java -jar /biosoftware/picard-tools-1.134/picard.jar CreateSequenceDictionary R=data/assembly.fasta O=data/assembly.dict

samtools faidx data/assembly.fasta

#过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates

java -jar /biosoftware/picard-tools-1.134/picard.jar MarkDuplicates I=./${i}/${i}_picard.bam O=./${i}/${i}.dedup.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=./${i}/${i}.metrics

#把所有的样本的bam文件都用HaplotypeCaller的GVCF模式分别分析变异

/data2/yanxu2016/gatk-4.0.0.0/gatk  --java-options -Xmx60g HaplotypeCaller  --emit-ref-confidence GVCF -R ./data/assembly.fasta -I ./${i}/${i}.dedup.bam -O ./gatk_out/${i}.HC.g.vcf

done

#合并gvcf

/data2/yanxu2016/gatk-4.0.0.0/gatk  --java-options -Xmx60g CombineGVCFs  -V ./gatk_out/fuben.HC.g.vcf -V ./gatk_out/muben.HC.g.vcf -V ./gatk_out/d1.HC.g.vcf -V ./gatk_out/s1.HC.g.vcf -R ./data/assembly.fasta -O ./gatk_out/Reseq_all_fmds.HC.vcf

#进行分型

/data2/yanxu2016/gatk-4.0.0.0/gatk  --java-options -Xmx60g GenotypeGVCFs -V Reseq_all_dfsm.HC.vcf -R ../data/assembly.fasta -O all_Reseq_genetype.HC.vcf

#使用samtools call variants 综合GATK和samtools的结果获得可靠变异

#之前安装的是1.1版本的samtools，1.1版本缺少对基因型信息统计的部分信息的筛选功能，所以重新装1.9版本，总是报错，查到github的project页面发现有些库不全。

https://github.com/samtools/samtools/blob/develop/INSTALL

./samtools-1.9/samtools mpileup -t AD,ADF,ADR,DP,SP -ugf ../data/assembly.fasta ../fuben/fuben.dedup.bam ../muben/muben.dedup.bam ../d1/d1.dedup.bam ../s1/s1.dedup.bam | bcftools call -vm > variants.samtools.raw1.vcf

未完待续