I have a CS background, but am a bioinformatics neophyte

I did a full genome sequence which provided me with ~100GB of files (SNP VCF, Indel VCF, BAM, Indel TBI, SNP TBI, BAM BAI, CNV VCF, CNV TBI, CV VCF, SV TBI, FASTQ R1, and FASTQ R2).

I'd like to search the genome for a pathogenic variant of a known gene. What tools do I need to do this, or where do I start?

  • 1
    $\begingroup$ Hi and welcome @ensnare Is only the gene known? or is the specific variant known as well? What format is the data you are searching for? $\endgroup$
    – Bioathlete
    Sep 28 '21 at 19:59
  • 1
    $\begingroup$ Well, does the vcf have gene annotation in it? or rsIDs? $\endgroup$
    – swbarnes2
    Sep 28 '21 at 23:10

I'm going to try and keep this answer from going overboard while still helping you find your feet. The good news is that it sounds like you have started off with some of this work already done for you, and based upon the files you have above, I will try to indicate what you have already and also where you might be able to go with this. The term you should be entering into Google as the generally accepted term for this pipeline is variant calling or possibly variant curation for the downstream steps.

  1. This process starts off with one or more files containing raw observed sequences from a sequencing device (often, but not always, Illumina) This files will usually be in the FASTQ format, but historically I have seen another format called QSEQ. You will likely never see this format and I feel old just bringing it up, but if you ever encounter it, just convert your files to FASTQ immediately. (Files entering here: FASTQ, possibly QSEQ if it was generated well before Carly Rae Jepson suggested that we call her, maybe)
  2. Although not a hard requirement, I recommend running something like fastqc or multiqc or some other form of quality control on your reads and making sure that the reads are of satisfactory quality with no anomalous metrics on them. You may or may not decide that you need to do some clean-up on the reads if you determine that there might be adapter contamination or that some reads are partially/entirely garbage quality. (Files entering here: FASTQ, reference genome in FASTA. Files generated: a quality metrics report and possibly a cleaned up FASTQ)
  3. Once you have your raw reads and consider them to be of satisfactory quality, they will need to be aligned to a reference genome. Read alignment is a critical, computationally intensive process that is required for many bioinformatics pipelines, meaning that there are many aligners out there that have been optimized for many different situations. If you are running a variant calling pipeline on Illumina short reads from genomic DNA, a current version of BWA should be a solid choice for you here. (Files entering here: FASTQ. Files generated: SAM and/or BAM files with aligned reads)
  4. Now that alignment is complete, we can start getting reads for the real variant calling process by doing any preprocessing needed. The most popular tool(s) for this historically were Picard Tools and GATK, but those two tool sets were used together for so long that they were just begging to be merged. The recent versions of GATK have made it so and Picard tools are now rolled into and called directly from GATK. If you are using an older GATK version, at a minimum you will want to recalibrate the base quality scores to try correcting biases that your sequencer had when assigning confidence scores to individual base calls (they all have this, trust me) and realigning reads around detected insertion/deletion (indel) sites. Newer versions of GATK don't need the indel realignment steps due to the haplotype calling algorithm going even further and reassembling reads around these complex variants. Those are just the minimal preprocessing steps, you may also need to assign read groups if they aren't already in the BAM files and possibly some other stuff (most likely GATK will let you know what's missing as you progress). (Files in: BAM, reference genome and "gold standard" variants in VCF format. Files out: preprocessed BAM)
  5. Actual variant calling can happen now, and this is likely going to happen using HaplotypeCaller (mutect2, if you're dealing with cancer/somatic changes). If this is your intention, you can look up the latest and greatest practices at https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels- for the most current guidance. (Files in: preprocessed BAM, reference genome. Files out: variant calls in VCF, possibly an active sites BAM file to visualize reassemblies around active indels)
  6. Filtering of variants should now happen. At this stage, you may be very application dependent, but essentially this process should assign some variants to being high confidence and others to being potentially technical artifacts. I have seen people skip this step either due to lack of interest or concern that it will filter out a real variant of interest, but I don't recommend that approach and have seen more people missing real variants due to being inundated with low confidence variants that should have been otherwise filtered out.
  7. Variant curation can now happen where one takes their high confidence variants and figures out (or usually has a computer program figure this out) which genes (if any) are being affected by the variant, what the chance will be to the protein coding as well as if the variant might affect splicing or some aspect of gene regulation or structure (such as an out of frame indel). It is often also useful here to compare to some database of mutation frequencies to filter out common SNP variants as potential cause for rare conditions.
  8. Identification of candidate mutations based upon the curation results. Identifying the mutations that are in genes that could plausibly control the observed phenotype, are uncommon, and are potentially causing pathogenic changes to the gene or its product (protein).

It looks like you already have VCF files, so at least the initial phase of variant calling has been finished for you already. I would recommend checking out something like Annovar or Varient Effect Predictor (VEP) that can run over your (hopefully filtered) VCF and assign a predicted gene and effect to a variant. Check both the SNP and Indel VCF files, as indels, while more rare, can have profound effects on a gene.

From your description, it sounds like you wanted to check a few genes in particular, and after the variant effect preduction, you should have the needed information to do this.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.