Determine whether a gene is totally missing

Question

[Canh]

I'm just a newbie in bioinformatics and need to do some testing on human genome.

After using some tools for SNP calling and annotation (platypus, snpeff), I can see a list of SNPs and which genes these SNPs affects, so I know these genes are modified compared to the reference genome. But how could I know if a gene is totally missing? Which files should I analyze? sam/bam/vcf... Or there's any tool doing this out there?

Thanks in advance.

Answer 1

[Biogeek]

More info needed...

Have you got a transcriptome? If so, best aligning to genome and doing gene counts. I'd recommend RSEM for counts . Gene which is missing in your reads will have a zero expression value at that locus.

[Canh]

I have one reference fasta file, 2 fastq files containing reads. I use bwa and samtools to align reads, then use platypus to call SNPs, then snpeff to annotate. That's all I did. So I need a transcriptome? Actually, I don't even know what a transcriptome is.

[WouterDeCoster]

It's not transcriptome sequencing if I understood correctly.

Answer 2

[DG]

Look at your alignment (SAM or BAM) that will show the actual sequencing reads mapped to the reference genome. A tool like IGV or IGB will let you view this and browse yo any area of interest (such as a gene) and visually see the reads as they are mapped. Are you dealing with Whole Genome or Exome sequencing data? In the case of Exome, depending on the gene, lack of data isn't conclusive for the gene actually being missing. Some genes and some exons in particular do not have very effective capture and are sometimes not sequenced, or have very little coverage.

[Canh]

I wanna write a script to do this, say input a gene (name, location) and output whether it is missing, modified, or the same as in the reference genome. Visual tools would not help me much. Are there any command line tools for this? I'm dealing with Exome, using Python. Thanks.

[igor]

You can use bedtools coverage to get the depth of coverage for specified regions (can be genes, exons, etc.): http://bedtools.readthedocs.io/en/latest/content/tools/coverage.html

[WouterDeCoster]

It's always bedtools saving the day

https://thescienceweb.wordpress.com/2015/03/25/most-bioinformaticians-to-be-replaced-by-bedtools/

[DG]

Yup, was going to be my suggestion as well

[WouterDeCoster]

What does it mean for a gene to be missing or modified? Missing = not covered by the sequencing? Modified = contains at least one variant? Does it have to be a coding variant?

I think you need think what exactly you want to achieve.

[DG]

This is important as well, what you mean by modified is important.

[DG]

See below for some caveats and suggestions. I'd also stress, again that absence of evidence is not the same thing as evidence of absence when it comes to exome sequencing. You get all sorts of exons and even whole genes with little to know coverage due to failure to capture or sequence, not necessarily because there is anything wrong with the actual gene in your sample.

What exactly are you looking for? Are you trying to identify copy-number variations, whole gene deletions, partial gene deletions, re-arrangements or other Structural Variants? What do you mean by modified? SNPs and small Indels? Or more like what I described above?