1.22-wes實戰-總結

annovar 註釋

mkdir annovar

cd annovar

#正確的檔案目錄：  

#/home/vip2/7E5241.chr1_2_raw.vcf

 #/home/qmcui/software/ANNOVAR/annovar/convert2annovar.pl

mkdir annovar && cd annovar

sample="7E5241"

/home/qmcui/software/ANNOVAR/annovar/convert2annovar.pl -format vcf4old /home/vip2/7E5241.chr1_2_raw.vcf >${sample}.annovar

第一小步先建立好annovar檔案

sample="7E5241"

/home/qmcui/software/ANNOVAR/annovar/annotate_variation.pl -buildver hg38 --outfile ${sample}.anno 7E5241.annovar /home/qmcui/software/ANNOVAR/annovar/humandb

？？？之後怎樣進行註釋怎樣檢視註釋結果就不會了？？？

IGV視覺化

bam建立索引

ls *fixed.bam|xargs -i samtools index {} 

？？？後面IGV視覺化就聽不懂了？？？

找變異

生成gvcf檔案

**這步很慢很慢，一定記得掛後臺，為了等它已經耽誤太多事兒了**

**nohup  命令 & 這樣就掛後臺了**

cat > sample.ID

7E5239

7E5240

7E5241

*把這三個數存到sample.ID中，做迴圈用*

ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"

snp="/home/qmcui/dbsnp_146.hg38.vcf.gz"

cat sample.ID | while read sample;do

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller -ERC GVCF -R $ref -I /home/qmcui/bqsr.bam/${sample}_bqsr.bam --dbsnp $snp -O ${sample}_raw.gvcf 1>${sample}.gvcf.log;done

經過漫長的等待會生成：7E5239_bqsr.bam_raw.gvcf 7E5240_bqsr.bam_raw.gvcf 7E5241_bqsr.bam_raw.gvcf

**並沒有報錯，但是一個半小時後生成的log日誌總是0k？？？怎麼回事？？？

合併gvcf，生成按染色體來找變異

seq 22 | while read id;do echo chr${id};done > bed.txt 

#把chr1一直到chr22存放在bed.txt檔案中（不會用for迴圈那就用while吧）

cat >> bed.txt 

chrM

chrX

chrY

#這個命令可以往bed.txt檔案中新增內容，我們把chrM、chrX、chrY新增到檔案中

ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"

snp="/home/qmcui/dbsnp_146.hg38.vcf.gz"

cat bed.txt | while read bed; do gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" -L $bed -R $ref -V 7E5239_bqsr.bam_raw.gvcf -V 7E5240_bqsr.bam_raw.gvcf -V 7E5241_bqsr.bam_raw.gvcf --genomicsdb-workspace-path gvcfs_${bed}.db

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" GenotypeGVCFs \ -R $ref -V gendb://gvcfs_${bed}.db --dbsnp $snp -O final_${bed}.vcf;done