Pythonバイオ?　Pythonバイオ/ツール?
3773　　　2019-06-19 (水) 16:57:45

Variant Calling 教科書向けサンプル　DRR006760 †

参照

↑

サンプルデータのダウンロード †

元研究は、Shimazaki H., Honda J, Naoi T., Namekawa M., Nakano I., Yazaki M., Nakamura K., Yoshida K., Ikeda S., Ishiura H., Fukuda Y., Takahashi Y., Goto J., Tsuji S. and Takiyama Y.: Autosomal-recessive complicated spastic paraplegia with a novel lysosomal trafficking regulator gene mutation. J Neurol Neurosurg Psychiatry. 2014 Sep;85(9):1024-8

DRY解析教本　p161 で使うDRR006760をダウンロード

アクセスは、登録データの検索　Web上で

http://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA000961

表示された右側のRUNからSRAをクリックして

/ddbj_database/dra/sralite/ByExp/litesra/DRX/DRX005/DRX005949/DRR006760/DRR006760.sra

をダウンロード。　　FASTQをクリックして２つのFASTQファイルをダウンロードしてもよさそうだ。

ダウンロードしたDRR006760.sraファイルを、dump-fastqで解凍分解。--split-filesを忘れずに。

dump-fastq --split-files DRR006760.sra

結果は

DRR006760_1.fastq
DRR006760_2.fastq

後で使うために、それぞれを圧縮。gzipでもbgzipでもよさそう。

gzip -c DRR006760_1.fastq
gzip -c DRR006760_2.fastq
↑

レファレンスゲノムのダウンロード †

DRY解析教本　p151 に従って、broadinstituteからhuman_g1k_v37_decoyをダウンロードする

wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/human_g1k_v37_decoy.dict.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/human_g1k_v37_decoy.fasta.fai.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/human_g1k_v37_decoy.fasta.gz
↑

インターバルリスト、既知のバリアントデータのダウンロード †

インターバルリスト

wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Broad.human.exome.b37.interval_list.gz

既知のバリアントデータ

wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/dbsnp_138.b37.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/dbsnp_138.b37.vcf.idx.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.idx.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.idx.gz
↑

クオリティチェック †

fastqc --nogroup DRR006760_1.fastq.gz
fastqc --nogroup DRR006760_2.fastq.gz
↑

トリミング †

java -jar -Xmx512m /usr/local/Trimmomatic/trimmomatic-0.38.jar \
  PE                     \
  -threads 32         \
  -phred33               \
  -trimlog log_DRR006760.txt \
  DRR006760_1.fastq.gz  \
  DRR006760_2.fastq.gz  \
  paired_DRR006760_1.trim.fastq.gz   \
  unpaired_DRR006760_1.trim.fastq.gz \
  paired_DRR006760_2.trim.fastq.gz   \
  unpaired_DRR006760_2.trim.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 \
  LEADING:20 \
  TRAILING:20 \
  SLIDINGWINDOW:4:15 \
  MINLEN:50

MINLENを50にしている（DRY解析での値を利用）。ILLUMINACLIPの妥当性は不明。DRY解析では指定していない。Illumina Adapter Sequences Document、　Illumina Adapter Sequences

↑

アラインメント †

レファレンスゲノム　human_g1k_v37_decoyにアラインメントする。マニュアル:ttps://ccb.jhu.edu/software/hisat2/manual.shtml

↑

HISAT-2でアラインメントする場合 †

まずhisat2-buildでレファレンスのインデックスを作る

nohup hisat2-build human_g1k_v37_decoy.fasta.gz human_g1k_v37_decoy-build.fasta &> nohup-hisat2-build.out

リファレンスの入力を圧縮済みgzで与えたところ、次のエラーで作成できず。

Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
Reference file does not seem to be a FASTA file
  Time to join reference sequences: 00:00:01
Total time for call to driver() for forward index: 00:50:11
Error: Encountered internal HISAT2 exception (#1)

リファレンスの入力を圧縮前のファイルにすると動作。

nohup hisat2-build human_g1k_v37_decoy.fasta human_g1k_v37_decoy-build.fasta &>nohup-hisat2-build-2.out &

マッピングは、

nohup hisat2 -p 32 --dta-cufflinks -x human_g1k_v37_decoy-build.fasta -1 paired_DRR006760_1.trim.fastq.gz -2 paired_DRR006760_2.trim.fastq.gz -S DRR006760_hisat2.sam >& nohup-hisat2-align.out

ソート・インデックス付け Samtools document

samtools sort -@16 -m4G DRR006760_hisat2.sam > DRR006760_hisat2_sorted.bam
samtools index DRR006760_hisat2_sorted.bam  DDR006760_hisat2_sorted.bai

GATK/picardによるリードペア情報の修復 GATK document

gatk FixMateInformation -I DRR006760_hisat2_sorted.bam -O DRR006760_hisat2_fixmate_sorted.bam -SO coordinate -AS true --CREATE_INDEX

重複リードの除去。GATK/picardで行う

gatk MarkDuplicates -I DRR006760_hisat2_fixmate_sorted.bam -M DRR006760_hisat2_fixmate.duplicated.metrics -O DRR006760_hisat2_fixmate_dedup_sorted.bam --ASSUME_SORTED TRUE --REMOVE_DUPLICATES TRUE
samtools index DRR006760_hisat2_fixmate_dedup_sorted.bam DRR006760_hisat2_fixmate_dedup_sorted.bai

これから、bcftoolsでmpileupし、VFCファイルを作る

bcftools mpileup -Ou -f human_g1k_v37_decoy.fasta DRR006760_hisat2_fixmate_dedup_sorted.bam | bcftools call -Ou -mv -o DRR006760_hisat2.call.vcf
↑

DRY教本通りBWAでアラインメントする場合 †

DRY教本通りのbwaによるアラインメントをやってみる。

bwa index -p human_g1k_v37_decoy human_g1k_v37_decoy.fasta.gz

bwa mem -t16 -M \
 -R "@RG\tID:FLOWCELLID\tSM:DDR006760\tPL:illumina\tLB:DRR006760_library_1" \
 human_g1k_v37_decoy.fasta \
 paired_DRR006760_1.trim.fastq.gz paired_DRR006760_2.trim.fastq.gz | \
 samtools view -@4 -bS - > DRR006760_aligned_reads.bam

結果はDRR006760_aligned_reads.bamに得られている。これをsortしindexしておく。

samtools sort -@16 -m4G DRR006760_aligned_reads.bam > DRR006760_aligned_reads_sorted.bam
samtools index DRR006760_aligned_reads_sorted.bam  DDR006760_aligned_reads_sorted.bai

GATK/picardによるリードペア情報の修復

gatk FixMateInformation -I DRR006760_aligned_reads_sorted.bam -O DRR006760_aligned_reads_fixmate_sorted.bam -SO coordinate -AS true --CREATE_INDEX

更にGATK/picardによる重複リードの除去

gatk MarkDuplicates -I DRR006760_aligned_reads_fixmate_sorted.bam -M DRR006760_aligned_reads_fixmate_duplicate.metrics -O DRR006760_aligned_reads_fixmate_dedup_sorted.bam --ASSUME_SORTED TRUE --REMOVE_DUPLICATES TRUE
samtools index DRR006760_aligned_reads_fixmate_dedup_sorted.bam DRR006760_aligned_reads_fixmate_dedup_sorted.bai

DRY解析教本ではこの次に「ローカルリアラインメント」でGATKのRealignerTargetCreatorとIndelRealignerが必要としているが、GATK4では外されているうえ、FORUM FORUMでは「要らないのでBest Practiceから外した。Just skip」と書いてある。

これから、bcftoolsでmpileupし、VFCファイルを作る

bcftools mpileup -Ou -f human_g1k_v37_decoy.fasta DRR006760_aligned_reads_fixmate_dedup_sorted.bam | bcftools call -Ou -mv -o DRR006760_aligned_reads.call.vcf

この状態ではBGZF （ブロックレベルで圧縮）してあるので、目で読みたければ伸長する必要アリ。

bctools view -o DRR006760_aligned_reads.call.readable.vcf DRR006760_aligned_reads.call.vcf

シーケンスクオリティの再評価　BaseRecalibrator

knownSitesとしてdbsnp_138.b37.vcfとMills_and_1000G_gold_standard.indels.b37.vcfを使って、再評価する。human_g1k_v37_decoy.dict、dbsnp_138.b37.vcf、Mills_and_1000G_gold_standard_indels.b37.vcfは解凍しておく必要があるらしい。

gatk BaseRecalibrator -I DRR006760_aligned_reads_fixmate_dedup_sorted.bam -R human_g1k_v37_decoy.fasta  --known-sites dbsnp_138.b37.vcf --known-sites Mills_and_1000G_Gold_standard.indels.b37.vcf -O DRR006760_recal.table

次のPrintReadsは、GATK4ではApplyBQSRに置き換えられているらしい。FORUM

gatk ApplyBQSR -R human_g1k_v37_decoy.fasta -I DRR006760_aligned_reads_fixmate_dedup_sorted.bam --bqsr-recal-file DRR006760_recal.table -O DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bam
samtools index DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bam DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bai

いよいよHaplotypeCallerでバリアント検出し、GVCFファイルを作る　（DRY解析入門p172）
準備で、Broad.human.exome.b37.interval.list.gz をungzipしておく。

gunzip -c Broad.human.exome.b37.interval.list.gz > Broad.human.exome.b37.interval.list
gatk --java-options "-Xmx4g" HaplotypeCaller -R human_g1k_v37_decoy.fasta -I DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bam -L Broad.human.exome.b37.interval.list --dbsnp dbsnp_138.b37.vcf -ERC GVCF -O DRR006760_raw_variants.g.vcf.gz -bamout bamout.bam

複数サンプルの集約処理　（ここでは１種類だけだが、と書いてある）

gatk --java-options "-Xmx4g" GenotypeGVCFs -R human_g1k_v37_decoy.fasta -V DRR006760_raw_variants.g.vcf.gz -O DRR006760_combined_genotyped.vcf.gz

フィルタリング処理のために、SelectVariantsを用いてSNV部分だけを抽出する

gatk SelectVariants -R human_g1k_v37_decoy.fasta -V DRR006760_combined_genotyped.vcf.gz --select-type-to-include SNP -O DRR006760_combined_genotyped_raw_snps.vcf

これは、pyVCFで簡単にできる

import vcf
infile = 'DRR006760_combined_genotyped.vcf'
outvcffile = 'DRR006760_combined_genotyped.snpsimple.vcf'

vcf_reader = vcf.Reader(open(infile, 'r'))
vcf_writer = vcf.Writer(open(outvcffile, 'w'), vcf_fh);
for i, snp_record in enumerate(vcf_reader):
    if snp_record.is_snp:
        vcf_writer.write_record(snp_record)
vcf_writer.close()
print('complete')

このis_snpの内容は、vcf.model._Recordのソースコードによると

   def is_snp(self):
       """ Return whether or not the variant is a SNP """
       if len(self.REF) > 1:
           return False
       for alt in self.ALT:
           if alt is None or alt.type != "SNV":
               return False
           if alt not in ['A', 'C', 'G', 'T', 'N', '*']:
               return False
       return True

同様に、is_indelやis_sv (structureal variant)、var_type（variant type）、var_subtypeなどがある。

SNV部分のフィルタリングを行う。GATKのフィルタリングは、FILTERフィールドにマークを付けるVariantFilteringと、実際にデータ中からレコードを削除するSelectVariantsの２段階で行える。（そうでなく、直接SelectVariantsで細かい条件を書いてしまうこともできそうだ）ここではDRY解析教本通りにやってみる。

gatk VariantFiltration -R human_g1k_v37_decoy.fasta \
  -V DRR006760_combined_genotyped_raw_snps.vcf \
  -O DRR006760_combined_genotyped_filtered_snps.vcf \
  --cluster-size 3 --cluster-window-size 10 \
  --filter-name "LowQD" --filter-expression "QD < 2.0" \
  --filter-name "HighFisherStrand" --filter-expression "FS > 60.0" \
  --filter-name "HighHaplotypeScore" --filter-expression "HaplotypeScore > 13.0" \
  --filter-name "LowRMSMappingQuality" --filter-expression "MQ < 40.0" \
  --filter-name "LowMQRankSum"  --filter-expression "MQRankSum < -12.5" \
  --filter-name "LowReadPosRankSum" --filter-expression "ReadPosRankSum < -8.0" \
  --filter-name "HardToValidate" \
     --filter-expression "MQ0 > 4 && ((MQ0 / (1.0 * DP)) > 0.1)" \
  --filter-name "VeryLowQual" --filter-expression "Qual < 30.0" \
  --filter-name "LowQual" --filter-expression "QUAL >= 30.0 && QUAL < 50.0"

INDEL部分だけのフィルタリング処理をするため、INDEL部分を取出す

gatk SelectVariants -R human_g1k_v37_decoy.fasta -V DRR006760_combined_genotyped.vcf.gz --select-type-to-include INDEL -O DRR006760_combined_genotyped_raw_indels.vcf.gz

上と同様に、VariantFilterationでクオリティの低いバリアントにマーク付けをする。

gatk VariantFiltration -R human_g1k_v37_decoy.fasta \
  -V DRR006760_combined_genotyped_raw_indels.vcf.gz \
  -O DRR006760_combined_genotyped_filtered_indels.vcf.gz \
  --filter-name "LowQD" --filter-expression "QD < 2.0" \
  --filter-name "HighFisherStrand" --filter-expression "FS > 200.0" \
  --filter-name "LowReadPosRankSum" --filter-expression "ReadPosRankSum < -20.0"\

CombineVariantsを使って、SNVのフィルタ結果とINDELのフィルタ結果を結合する。CombineVariantsはgatkで呼べるサブコマンドに入っていないので、直接javaから起動する。

java -Xmx4g -jar /usr/local/src/gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar \
  -T CombineVariants -R human_g1k_v37_decoy.fasta  \
  --assumeIdenticalSamples \
  -V:SNP DRR006760_combined_genotyped_filtered_snps.vcf.gz \
  -V:INDEL DRR006760_combined_genotyped_filtered_indels.vcf.gz \
  -oDRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf.gz \
  -genotypeMergeOptions UNIQUIFY

上記はうまく動かないようだ。代わりに以下を使う。

gatk MergeVcfs -I DRR006760_combined_genotyped_filtered_snps.vcf \
  -I DRR006760_combined_genotyped_filtered_indels.vcf \
  -O DRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf.gz

（gzは処理できないのでungzipしておく）

SelectVariantsでマークしたデータを除去する。

gatk SelectVariants -R human_g1k_v37_decoy.fasta \
  -V DRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf.gz \
  --exclude-filtered --exclude-non-variants \
  -O DRR006760_combined_genotyped_filtered_snps_indels_mixed.PASS.vcf

p182 バリアントのアノテーション Annovarを使っている。

代わりに、SpnEffを使ってみよう。Javaベース。マニュアル

java -Xmx4g -jar /usr/local/snpEff.jar GRCh37.75 DRR006760_combined_genotyped_filtered_snps_indels_mixed.PASS.vcf > DRR006760_PASS.ann.vcf

生成されたDRR006760_PASS.ann.vcfの先頭部分

4       190878597       .       C       T       1554.6  PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=-4.162e+00;DP=235;ExcessHet=3.0103;FS=120.742;MLEAC=1;MLEAF=0.500;MQ=43.66;MQRankSum=-1.121e+01;QD=6.64;ReadPosRankSum=3.14;SOR=5.283;ANN=T|synonymous_variant|LOW|FRG1|ENSG00000109536|transcript|ENST00000226798|protein_coding|6/9|c.477C>T|p.Cys159Cys|699/1074|477/777|159/258||,T|synonymous_variant|LOW|FRG1|ENSG00000109536|transcript|ENST00000524583|protein_coding|5/7|c.93C>T|p.Cys31Cys|540/755|93/308|31/101||WARNING_TRANSCRIPT_INCOMPLETE,T|synonymous_variant|LOW|FRG1|ENSG00000109536|transcript|ENST00000531991|protein_coding|7/7|c.288C>T|p.Cys96Cys|678/738|288/348|96/115||WARNING_TRANSCRIPT_NO_STOP_CODON,T|upstream_gene_variant|MODIFIER|FRG1|ENSG00000109536|transcript|ENST00000507103|nonsense_mediated_decay||c.-3306C>T|||||3306|WARNING_TRANSCRIPT_NO_START_CODON,T|upstream_gene_variant|MODIFIER|FRG1|ENSG00000109536|transcript|ENST00000505327|processed_transcript||n.-3640C>T|||||3640|,T|downstream_gene_variant|MODIFIER|FRG1|ENSG00000109536|transcript|ENST00000533157|nonsense_mediated_decay||c.*15543C>T|||||2292|,T|non_coding_transcript_exon_variant|MODIFIER|FRG1|ENSG00000109536|transcript|ENST00000514482|processed_transcript|6/7|n.463C>T||||||  GT:AD:DP:GQ:PL  0/1:161,73:234:99:1562,0,5571
6       41774685        .       C       G       229.6   PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=-2.290e+00;DP=51;ExcessHet=3.0103;FS=69.864;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=4.59;ReadPosRankSum=-5.193e+00;SOR=5.992;ANN=G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000373009|protein_coding|1/5|c.37G>C|p.Ala13Pro|37/2067|37/2067|13/688||,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000394253|protein_coding|3/7|c.37G>C|p.Ala13Pro|367/9034|37/2067|13/688||,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000373010|protein_coding|6/10|c.37G>C|p.Ala13Pro|532/2446|37/1758|13/585||,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000373006|protein_coding|4/7|c.37G>C|p.Ala13Pro|259/3177|37/1923|13/640||,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000297229|protein_coding|2/5|c.37G>C|p.Ala13Pro|145/3062|37/1923|13/640||,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000437061|protein_coding|2/2|c.37G>C|p.Ala13Pro|477/553|37/113|13/36||WARNING_TRANSCRIPT_INCOMPLETE,G|missense_variant|MODERATE|USP49|ENSG00000164663|transcript|ENST00000423567|protein_coding|4/4|c.37G>C|p.Ala13Pro|744/752|37/45|13/14||WARNING_TRANSCRIPT_NO_STOP_CODON,G|upstream_gene_variant|MODIFIER|USP49|ENSG00000164663|transcript|ENST00000448078|nonsense_mediated_decay||c.-1203G>C|||||1202|WARNING_TRANSCRIPT_NO_START_CODON       GT:AD:DP:GQ:PL  0/1:31,19:50:99:237,0,831
6       71187020        .       A       C       260.6   PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=-4.040e+00;DP=88;ExcessHet=3.0103;FS=122.727;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=2.96;ReadPosRankSum=-4.470e+00;SOR=7.222;ANN=C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000418814|protein_coding|8/22|c.527A>C|p.His176Pro|1141/6415|527/4548|176/1515||,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000370479|protein_coding|7/22|c.398A>C|p.His133Pro|916/5676|398/3909|133/1302||,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000505769|protein_coding|7/22|c.527A>C|p.His176Pro|934/4806|527/3288|176/1095||,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000515323|protein_coding|7/14|c.527A>C|p.His176Pro|938/1926|527/1515|176/504||WARNING_TRANSCRIPT_NO_STOP_CODON,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000457062|protein_coding|6/21|c.398A>C|p.His133Pro|668/5282|398/3909|133/1302||,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000361499|protein_coding|8/22|c.527A>C|p.His176Pro|871/5042|527/3960|176/1319||,C|missense_variant|MODERATE|FAM135A|ENSG00000082269|transcript|ENST00000505868|protein_coding|5/18|c.527A>C|p.His176Pro|527/5173|527/4410|176/1469||,C|synonymous_variant|LOW|FAM135A|ENSG00000082269|transcript|ENST00000194672|nonsense_mediated_decay|6/21|c.393A>C|p.Thr131Thr|437/5198|393/444|131/147||,C|upstream_gene_variant|MODIFIER|FAM135A|ENSG00000082269|transcript|ENST00000425415|retained_intron||n.-4816A>C|||||4816|,C|non_coding_transcript_exon_variant|MODIFIER|FAM135A|ENSG00000082269|transcript|ENST00000393299|retained_intron|5/12|n.484A>C||||||     GT:AD:DP:GQ:PL  0/1:52,36:88:99:268,0,1339

このデータをpyVCFで読んでみる。

import vcf
# SnpEffで得られた　DRR006760_PASS.ann.vcf　をpyVCFで読んでみる
infile = 'DRR006760_PASS.ann.vcf'
outfile = 'mydata.txt'
with open(outfile, 'w') as outfh:
    vcf_fh = vcf.Reader(open(infile, 'r'))
    for i, snp_record in enumerate(vcf_fh):
        if i>5:
            break
        pl_list = []
        print(dir(snp_record))
        print('INFO', snp_record.INFO)
        print('ANN', snp_record.INFO['ANN'])
        for sample in snp_record.samples:
            #print(dir(sample.data))
            pl_list.append(str(sample['PL']))
        print(pl_list)
            
        new_record = [str(snp_record.CHROM),
                        str(snp_record.POS),
                        str(snp_record.ID),
                        str(snp_record.REF),
                        str(snp_record.ALT)]
            
        new_record.extend(pl_list)
            
        #outfh.write('\t'.join(new_record) + '\n')
print('complete')

また、ファイル　snpEff_summary.html　と　snpEff_genes.txt　が得られた。

snpEff_summary.htmlはWeb表示になっていて、

SnpEff_html1.png
SnpEff_html2.png
SnpEff_html3.png
SnpEff_html4.png
SnpEff_html5.png
SnpEff_html6.png
SnpEff_html7.png
SnpEff_html8.png
SnpEff_html9.png
SnpEff_html10.png
SnpEff_html10-2.png
SnpEff_html11.png

snpEff_genes.txt はタブ区切りで上図に対応するデータが変異ごとに表にされている。

SnpEff_genes2.png
↑

VCFのフィルタリングを自分バージョンで †

import vcf

# pyVCFでSNP抽出し、フィルタを掛けてみる。pyVCFのテンプレートにフィルタを足す方法もあるが、ここは自前とする。
# DRY解析教本に出ているフィルタ（以下の通り）と同じことをする
# gatk VariantFiltration -V DRR006760_combined_genotyped_raw_snps.vcf -O DRR006760_combined_genotyped_filtered_snps.vcf \
#   --cluster-size 3 --cluster-window-size 10 \
#   --filter-name "LowQD" --filter-expression "QD < 2.0" \  QDはQuolity Score(6列目)/NumberOfReads
#   --filter-name "HighFisherStrand" --filter-expression "FS > 60.0" \
#   --filter-name "HighHaplotypeScore" --filter-expression "HaplotypeScore > 13.0" \  
#   --filter-name "LowRMSMappingQuality" --filter-expression "MQ < 40.0" \
#   --filter-name "LowMQRankSum" --filter-expression "MQRankSum < -12.5" \
#   --filter-name "LowReadPosRankSum" --filterExpression "ReadPosRankSum < -8.0" \
#   --filter-name "HardToValidate" --filter-expression "MQ0 > 4 && ((MQ0 / (1.0 * DP)) > 0.1)" \
#   --filter-name "VeryLowQual" --filter-expression "Qual < 30.0" \
#   --filter-name "LowQual" --filter-expression "QUAL >= 30.0 && QUAL < 50.0" 
# GATKのページhttps://software.broadinstitute.org/gatk/documentation/article?id=6925では
#    QD, FS, MQ, MQRankSum and ReadPosRankSum.を使うと書いている。
def filtered1(rec):
# もしフィルタ条件を満たせばパス（True）、満たさなければドロップ（False）
    
   if 'QD' in rec.INFO.keys():
       QD = (rec.INFO['QD'] < 2.0)
   else:
       QD = False  # QD記述がなければfilterしない
   if 'FS' in rec.INFO.keys():
       FS = (rec.INFO['FS'] > 60.0)
   else:
       FS = False
   if 'HighHaplotypeScore' in rec.INFO.keys():
       HighHaplotypeScore = (rec.INFO['HighHaplotypeScore'] > 60.0)
   else:
       HighHaplotypeScore = False
   if 'MD' in rec.INFO.keys():
       MQ = (rec.INFO['MQ'] < 40.0)
   else:
       MQ = False
   if 'MQRankSum' in rec.INFO.keys():
       MQRankSum = (rec.INFO['MQRankSum'] < -12.0)
   else:
       MQRankSum = False
   if 'ReadPosRankSum' in rec.INFO.keys():
       ReadPosRankSum = (rec.INFO['ReadPosRankSum'] < -8.0)
   else:
       ReadPosRankSum = False
   LowQual = (rec.QUAL < 50)

   filter = QD or FS or HighHaplotypeScore or MQ or MQRankSum or ReadPosRankSum or LowQual
   return(filter)

def filtered2(rec):
# もしフィルタ条件を満たせばパス（True）、満たさなければドロップ（False）
    
    if 'QD' in rec.INFO.keys():
        QD = (rec.INFO['QD'] < 2.0)
    else:
        QD = False  # QD記述がなければfilterしない
    if 'FS' in rec.INFO.keys():
        FS = (rec.INFO['FS'] > 200.0)
    else:
        FS = False
    if 'ReadPosRankSum' in rec.INFO.keys():
        ReadPosRankSum = (rec.INFO['ReadPosRankSum'] < -20.0)
    else:
        ReadPosRankSum = False
    filter = QD or FS or ReadPosRankSum
    return(filter)

infile = 'DRR006760_combined_genotyped.vcf'
#outfile = 'DRR006760_combined_genotyped.vcf.txt'
outvcffile = 'DRR006760_combined_genotyped.snpsimple.vcf'

vcf_fh = vcf.Reader(open(infile, 'r'))
vcfwriter = vcf.Writer(open(outvcffile, 'w'), vcf_fh);
for i, snp_record in enumerate(vcf_fh):
    #if (not snp_record.is_snp) or filtered(snp_record):
    #    print(i, snp_record)
    if (snp_record.is_snp and not filtered1(snp_record)) or (snp_record.is_indel and not filtered2(snp_record)):
        vcfwriter.write_record(snp_record)
vcfwriter.close()
#vcf_fh.close()
print('complete')

この出力を

java -Xmx4g -jar /usr/local/snpEff.jar GRCh37.75 DRR006760_combined_genotyped.snpsimple.vcf > DRR006760_snpsimple.ann.vcf

結果は、

wc -l snpEff_genes-0609.txt
64936 snpEff_genes-0609.txt
に対して
wc -l snpEff_genes-0609-snpsimple.txt
65735 snpEff_genes-0609-snpsimple.txt
SnpEff-snpsimple_html1.png
SnpEff-snpsimple_html7.png

GATKでのフィルタに比べると少し多く残っていることが分かる。 何かうまく行っていない。

↑

追加　2019-06-12 †

マッピング

NIGのマッピングツール比較 (2013/6におけるオプション比較)

バリアントコール・パイプラインのサーベイ/比較

論文： A Comparison of Variant Calling Pipelines Using Genome in a Bottle as a Reference
Hindawi Publishing Corporation BioMed Research International Volume 2015, Article ID 456479, 11 pages http://dx.doi.org/10.1155/2015/456479

論文： A beginners guide to SNP calling from high-throughput DNA-sequencing data
Altmann A, Weber P, Bader D, Preuss M, Binder EB, Müller-Myhsok B. Hum Genet. 2012 Oct;131(10):1541-54.

日本語要約　review article要約　SNPs callingビギナーズガイド - macでインフォマティクス

FreeBayes コーリングソフト

FreeBayes ekg/freebayes: Bayesian haplotype-based genetic polymorphism discovery and genotyping.
[1207.3907v2] Haplotype-based variant detection from short-read sequencing (2012 version)
Haplotype-based variant detection from short-read sequencing (2016 version, slightly updated)

freebayes インストール 簡単な使い方 – バイオインフォ 道場 [bioinfo-Dojo]

解説

ハプロタイプ解析における推定・ブロック同定および関連解析（冨田 誠）

Base Quality Score Recalibration (BQSR)
Unfortunately the scores produced by the machines are subject to various sources of systematic (non-random) technical error, leading to over- or under-estimated base quality scores in the data. Some of these errors are due to the physics or the chemistry of how the sequencing reaction works, and some are probably due to manufacturing flaws in the equipment.

VCFや他のアノテーションを結合してアノテーションを作る？

Vcfanno: fast, flexible annotation of genetic variants (2016)
brentp/vcfanno: annotate a VCF with other VCFs/BEDs/tabixed files

添付ファイル: fileSnpEff-snpsimple_html7.png 873件 [詳細] fileSnpEff-snpsimple_html1.png 847件 [詳細] fileSnpEff_html10-2.png 819件 [詳細] fileSnpEff_genes2.png 786件 [詳細] fileSnpEff_html11.png 788件 [詳細] fileSnpEff_html10.png 790件 [詳細] fileSnpEff_html9.png 816件 [詳細] fileSnpEff_html8.png 751件 [詳細] fileSnpEff_html7.png 808件 [詳細] fileSnpEff_html6.png 754件 [詳細] fileSnpEff_html5.png 741件 [詳細] fileSnpEff_html4.png 706件 [詳細] fileSnpEff_html3.png 736件 [詳細] fileSnpEff_html2.png 790件 [詳細] fileSnpEff_html1.png 780件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-06-19 (水) 16:57:45 (1379d)