#文件准备
ref=dog_chr5.fa
R1=BD225_TAGCTT_L007_R1_001.pe.fq.gz
R2=BD225_TAGCTT_L007_R1_002.pe.fq.gz
#工具准备
bwa
samtools
bamdst
GATK3

#step1
ln -s example/data/Canis_lupus_familiaris/BD225_TAGCTT_L007_R1_002.pe.fq.gz
ln -s example/data/Canis_lupus_familiaris/BD225_TAGCTT_L007_R1_001.pe.fq.gz
cp /home/lixin/example/data/Canis_lupus_familiaris/dog_chr5.fa.gz .
gzip -d dog_chr5.fa.gz
bwa index -a bwtsw dog_chr5.fa
/home/lixin/example/program/samtools-1.3.1/bin/samtools faidx dog_chr5.fa
/home/lixin/example/program/samtools-1.3.1/bin/samtools dict dog_chr5.fa -o dog_chr5.dict

#step5 get read group info
R1=BD225_TAGCTT_L007_R1_001.pe.fq.gz
SM=$(echo  | cut -d"_" -f1)
LB=$(echo $R1 | cut -d"_" -f1,2)                                   ##library ID
PL="Illumina"                                                           ##platform (e.g. illumina, solid)
RGID=$(zcat $R1 | head -n1 | sed 's/:/_/g' |cut -d "_" -f1,2,3,4)       ##read group identifier
PU=$RGID.$LB                                                            ##Platform Unit
echo -e "@RG\tID:$RGID\tSM:$SM\tPL:$PL\tLB:$LB\tPU:$PU"
#"@RG\tID:@D3VG1JS1_213_C7R8WACXX_7\tSM:BD225\tPL:Illumina\tLB:BD225_TAGCTT\tPU:@D3VG1JS1_213_C7R8WACXX_7.BD225_TAGCTT"

#step6 bwa mem and convert default sam output into bam and sorted bam
ref=dog_chr5.fa
R1=BD225_TAGCTT_L007_R1_001.pe.fq.gz
R2=BD225_TAGCTT_L007_R1_002.pe.fq.gz
bwa mem -t 4 -M -R "@RG\tID:@D3VG1JS1_213_C7R8WACXX_7\tSM:BD225\tPL:Illumina\tLB:BD225_TAGCTT\tPU:@D3VG1JS1_213_C7R8WACXX_7.BD225_TAGCTT" $ref BD225_TAGCTT_L007_R1_001.pe.fq.gz BD225_TAGCTT_L007_R1_002.pe.fq.gz | /home/lixin/example/program/samtools-1.3.1/bin/samtools view -bS - |/home/lixin/example/program/samtools-1.3.1/bin/samtools sort -l 9 -m 800M -@ 4 -T sorted -o L007.sorted.bam -

wget 'ftp://ftp.ensembl.org/pub/release-89/variation/vcf/canis_familiaris/Canis_familiaris.vcf.gz' -O canis_familiaris.vcf.gz
gzip -dc canis_familiaris.vcf.gz | grep -E "^#|^5" | sed 's/^5/chr5/' >canis_fam_chr5.vcf
ref=dog_chr5.fa
knownSites=canis_fam_chr5.vcf
sample=L007.sorted.bam
/home/lixin/example/program/samtools-1.3.1/bin/samtools index L007.sorted.bam
java -Xmx2g -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R $ref -I $sample -knownSites $knownSites -o L007.1st.table
java -Xmx2g -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R $ref -I $sample -knownSites $knownSites -BQSR L007.1st.table -o L007.2nd.table
java -Xmx2g -jar GenomeAnalysisTK.jar -T PrintReads -R $ref -I $sample -BQSR L007.2nd.table -o L007.recal.bam
java -Xmx2g -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R $ref --dbsnp $knownSites -I $sample --emitRefConfidence GVCF -nct 3 -o L007.g.vcf
java -Xmx2g -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R $ref --dbsnp $knownSites --variant L007.g.vcf -o raw_variants.vcf
wget https://github.com/shiquan/bamdst/archive/master.zip
unzip master.zip
cd bamdst-master
make
cd ../
./bamdst-master/bamdst -p chr5.bed -o ./coverage  L007.recal.bam

java -jar  GenomeAnalysisTK.jar  -R  $ref -T SelectVariants -nt 8 --variant raw_variants.vcf  -o snv.raw.vcf -selectType SNP
java -jar  GenomeAnalysisTK.jar   -R  $ref -T VariantFiltration --variant snv.raw.vcf -o snv.filter.vcf --filterExpression "QD<2.0 || MQ<40.0 || FS>60.0 " --filterName "StandardFilter"
java -jar  GenomeAnalysisTK.jar  -R  $ref -T SelectVariants -nt 8 --variant raw_variants.vcf  -o indel.raw.vcf -selectType INDEL
java -jar  GenomeAnalysisTK.jar   -R  $ref -T VariantFiltration --variant indel.raw.vcf -o indel.filter.vcf --filterExpression "QD<2.0 || FS>200.0 " --filterName "StandardFilter"
awk '{if($1~/^chr5/){print $1}}' snv.raw.vcf|wc -l
awk '{if($1~/^chr5/){print $1}}' indel.raw.vcf|wc -l
##平均覆盖深度：0.25 参考序列： 5.50%
##SNP个数：5411 INDEL个数：1395