Frage:
Keine Variante mit GATK 4.0 HaplotypeCaller gefunden
Ammar Sabir Cheema
2018-07-04 00:49:04 UTC
view on stackexchange narkive permalink

Ich mache einen Variantenaufruf für RNA-seq-Datensätze aus Weizen, der hexaploide ist. Die BAM-Dateien (Binary Alignment) wurden mit STAR Version 2.6.0c erstellt und der Variantenaufruf wurde mit GATK 4.0 HaplotypeCaller durchgeführt. Die gesamte Pipeline ist wie folgt :

Ausrichtung durch STAR 2.6.0c : 

  STAR --runMode alignReads --runThreadN 20 --outSAMtype BAM SortedByCoordinate --sjdbOverhang 100 --genomeDir Weizen --readFilesIn G1_cleaned_R1.fastq G1_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript Wheat.gff3

Dann wurde AddOrReplaceReadGroups von picard verwendet mit Standardeinstellungen: 

  AddOrReplaceReadGroups \ I = input.bam \ O = Read_groups_added.bam \ RGID = 4 \ RGLB = lib1 \ RGPL = illumina \ RGPU = unit1 \ RGSM = 20

Schließlich wurde GATK 4.0 haplotypeCaller angewendet: 

  gatk HaplotypeCaller --native-pair-hmm-threads 28 \ -I Read_groups_added.bam -R Weizen_ref.fa -O fina l.vcf

Der Befehl wurde wie folgt erfolgreich ausgeführt: 

  20: 33: 01.018 INFO ProgressMeter - TGACv1_scaffold_470693_5DS: 301 102.2 45108670 441320.020: 33: 11.039 INFO ProgressMeter - TGACv1_scaffold_724195_U: 301 102,4 45.111.270 440625,520: 33: 21,072 INFO ProgressMeter - TGACv1_scaffold_724579_U: 301 102,5 45.113.880 439932,420: 33: 31,108 INFO ProgressMeter - TGACv1_scaffold_248869_3B: 301 102,7 45.116.440 439240,920: 33: 41,143 INFO ProgressMeter - TGACv1_scaffold_248918_3B: 301 102,9 45.118.990 438551,620: 33 : 51.150 INFO ProgressMeter - TGACv1_scaffold_306137_4AL: 301 103.0 45121550 437866.720: 34: 01.176 INFO ProgressMeter - TGACv1_scaffold_193559_2DS: 301 103.2 45124100 437182.5
20: 34: 03,134 INFO HaplotypeCaller - 100.112.154 read (s) gefiltert, durch: ((((((((MappingQualityReadFilter UND MappingQualityAvailableReadFilter) UND MappedReadFilter) UND NotSecondaryAlignmentReadFilter) UND NotDuplicateReadFilter) UND PassesVendorQualityCheckReadFilter) UND NonZeroReferenceLengthAlignmentReadFilter) UND GoodCigarReadFilter) UND WellformedReadFilter) 100.112.154 read (s) gefiltert nach: ((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassedVendorQualityCheckReadFilter) MappingQualityReadFilter UND MappingQualityAvailableReadFilter) UND MappedReadFilter) UND NotSecondaryAlignmentReadFilter) UND NotDuplicateReadFilter) UND PassesVendorQualityCheckReadFilter) UND NonZeroReferenceLengthAlignmentReadFilter MappingQualityAvailableReadFilter ND) UND MappedReadFilter) UND NotSecondaryAlignmentReadFilter) UND NotDuplicateReadFilter) UND PassesVendorQualityCheckReadFilter) 100112154 Lese (n) gefiltert nach: ((((MappingQualityReadFilter UND MappingQualityAvailableReadFilter) UND MappedReadFilter) UND NotSecondaryAlignmentReadFilter) UND NotDuplicateReadFilter) 100112154 Lese (n) gefiltert nach: (( (MappingQualityReadFilter UND MappingQualityAvailableReadFilter) UND MappedReadFilter) UND NotSecondaryAlignmentReadFilter) 100112154 Lese (n) gefiltert nach: ((MappingQualityReadFilter UND MappingQualityAvailableReadFilter) UND MappedReadFilter) 100112154 Lese (n) gefiltert nach: (MappingQualityReadFilter UND MappingQualityAvailableReadFilter) 21066679 Lese (n) gefiltert nach: MappingQualityReadFilter 79045475 Lesevorgänge gefiltert nach: MappingQualityAvailableReadFilter 20: 34: 03.134 INFO ProgressMeter - TGACv1_scaffold_725922_U: 301 103.2 45124592 437049.1
20: 34: 03.134 INFO ProgressMeter - Durchlauf abgeschlossen. Verarbeitete 45124592 Gesamtregionen in 103,2 Minuten.20: 34: 04.749 INFO VectorLoglessPairHMM - Zeitaufwand für die Einrichtung des JNI-Aufrufs: 0.020: 34: 04.749 INFO PairHMM - Gesamtberechnungszeit in PairHMM computeLogLikelihoods (): 0.020: 34: 04.749 INFO SmithW Rechenzeit in Java Smith-Waterman: 0,00 Sek. 20: 34: 04,749 INFO HaplotypeCaller - Motor abstellen [3. Juli 2018 20:34:04 PKT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller erledigt. Verstrichene Zeit: 103,99 Minuten.Runtime.totalMemory () = 4789370880

Die generierte vcf-Datei, dh final.vcf, enthält jedoch keine der folgenden Varianten: 

  [dellemc @ localhost bilal] $ tail final.vcf ## Contig = <ID = TGACv1_scaffold_725916_U, Länge = 500> ## Contig = <ID = TGACv1_scaffold_725917_U, Länge = 500> ## Contig = <ID = TGACv1_scaffold_725918_U, Länge = 500> ## Contig = <ID = TGACv1_scaffold_725919_U, Länge = ## 500> Contig = <ID = TGACv1_scaffold_725920_U, Länge = ## 500> Contig = <ID = TGACv1_scaffold_725921_U, Länge = ## 500> Contig = <ID = TGACv1_scaffold_725922_U, Länge = ## 500> Contig = <ID = TGACv1_scaffold_731344_7BL, Länge = 500> ## source = HaplotypeCaller # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 20

Liegt ein Fehler im Befehl haplotypeCaller vor oder fehlt ein Argument?

Alle diese Rechenzeiten von 0,0 scheinen etwas verdächtig.
Ich habe die Frage bearbeitet, um weitere Details hinzuzufügen.
@swbarnes können Sie über die Art des Verdachts berichten?
Zwei antworten:
sjcockell
2018-07-04 15:55:54 UTC
view on stackexchange narkive permalink

Es fehlen Schritte des Workflows "Best Practices", die meiner Meinung nach dazu führen, dass die meisten Lesevorgänge vom Haplotyp-Aufrufer gefiltert werden. Beachten Sie die folgenden Zeilen in der Ausgabe: 

  21066679 Lesevorgänge gefiltert nach: MappingQualityReadFilter 79045475 Lesevorgänge gefiltert nach: MappingQualityAvailableReadFilter

Das sind> 100 Millionen Lesevorgänge, die vom HC gefiltert werden.

Speziell sind Sie Fehlender SplitNCigarReads , der Lesevorgänge trennt, die eine Spleißverbindung überbrücken (in der STAR-Ausrichtung mit Ns aufgefüllt), und die Zuordnungsqualität neu zuweist. Derzeit ist dieses Element des Workflows in GATK4 nicht validiert oder verfügbar (der Aufruf von RNA-Seq-Varianten ist in neueren Versionen von GATK, AFAIK weitgehend aufgegeben), aber in den Best Practices verwendet WDL GATK3.5 - siehe https: //github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl, Zeile 440-477.

Es fehlen auch andere Elemente des empfohlenen Workflows, aber ich denke, dies ist das kritische Bit. Ich werde hier nicht auf die Ungeeignetheit von RNA-Seq-Daten für Variantenaufrufe eingehen, da mir nur allzu bewusst ist, dass es oft nicht viele Optionen gibt.

Was ist mit der Option --limitOutSJcollapsed im STAR-Befehl gemeint, der in dem von Ihnen vorgeschlagenen Link angegeben ist? -variant-calling.wdl # L345
Eine "reduzierte" Junction wird während eines STAR-Laufs in die "SJ.out.tab" -Datei geschrieben (die Reduktion bezieht sich hier auf das Entfernen doppelter Junctions). Diese Option begrenzt lediglich die maximale Anzahl dieser Junctions, die geschrieben werden. Der Standardwert ist 1.000.000.
Matt Bashton
2018-07-13 20:04:11 UTC
view on stackexchange narkive permalink

Neben den von @sjcockell erwähnten Problemen haben Sie ein weiteres potenzielles Problem, bei dem Sie einen hexaploiden Organismus haben. Standardmäßig ist der HaplotypeCaller (HC) fest auf eine Ploidie von 2 verdrahtet. Sie müssen diesen Wert auf 6 setzen, um Ihren Anforderungen zu entsprechen Organismen Ploidie unter Verwendung von -Ploidie 6 als Argument. Da die Ploidie ein Schlüsselbegriff in den zugrunde liegenden Bayes'schen Statistiken ist, die der HC verwendet (eine Erklärung finden Sie unter diesem Link), ist es wahrscheinlich, dass einige der niedrigeren Allelfrequenzen, die durch Hexaploidie erzeugt werden, einige gewonnene Varianten bedeuten Es wird nicht genannt, da es für einen diploiden Organismus statistisch unwahrscheinlich wäre.

Vielleicht möchten Sie auch das GATK-Supportforum in Betracht ziehen, da es in GATK-Angelegenheiten sehr hilfreich ist: https://gatkforums.broadinstitute.org/gatk/categories/gatk-support-forum


ⓘ
Diese Fragen und Antworten wurden automatisch aus der englischen Sprache übersetzt.Der ursprüngliche Inhalt ist auf stackexchange verfügbar. Wir danken ihm für die cc by-sa 4.0-Lizenz, unter der er vertrieben wird.
Loading...