0
$\begingroup$

Greetings,

I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:

1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3).

2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).

3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15).

4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).

5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).

6- Calling variants again using BCFtools.

7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools.

8- Filtering variants using BCFtools.

9- Select variants SNPs and Indels in two separated files using GATK.

10- Predicting effects using SnpEff.

11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.

Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence?

In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not?

Your response is highly appreciated,

Kind regards,

Improve this question
$\endgroup$
3
  • $\begingroup$ One thing that could be highly beneficial to this perfectly reasonable pipeline is blasting against the human genome and seeing whether it is conserved, there are any gnomAD mutations or worse clinvar mutations —of particular interest are the pLI scores. We know a lot about humans, so why not take advantage of it. $\endgroup$
    – Matteo Ferla
    Commented Nov 13, 2019 at 13:44
  • $\begingroup$ Do you mean blasting variants vcf file to the human genome? can you elaborate more? $\endgroup$
    – Mohamed Moawad
    Commented Nov 14, 2019 at 8:36
  • $\begingroup$ Yes, it's just an idea of a potential annotation. gnomAD is a very large DB of human variants and uses Ensembl ids, so, yes, a local blastp the camel protein canonical transcripts to the Ensembl fasta list of canonical transcripts and then the creation of a mapping function to see if the VCF would be pathogenic were it human or in a site also variable in humans. $\endgroup$
    – Matteo Ferla
    Commented Nov 15, 2019 at 17:37

2 Answers 2

Reset to default
1
$\begingroup$

This is close to what many people do. May I add:

  • between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.

  • filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold. For that purpose, some time ago I made this and this.

Improve this answer
$\endgroup$
0
$\begingroup$

It looks like you mostly followed GATK? It looks reasonable, though a lot will depend on how good your genome and gtf file are.

Improve this answer
$\endgroup$
1
  • $\begingroup$ You seems right about the genome and the gtf file. I'm following a reference paper methods and the pipeline produced reasonable number of SNPs and Indels, however, when it comes to predict the effects of those variants SnpEff gives me different predictions from the reference paper. Note that the genome is a scaffold level and the original gff file doesn't contain intron and intergenic regions. $\endgroup$
    – Mohamed Moawad
    Commented Nov 14, 2019 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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