<p>This is close to what many people do. May I add:</p>

<ul>
<li><p>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 <a href="https://github.com/MatteoSchiavinato/Utilities/blob/master/add-read-groups.py" rel="nofollow noreferrer">this</a>. </p></li>
<li><p>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 <a href="https://github.com/MatteoSchiavinato/Utilities/blob/master/SNP-dogg" rel="nofollow noreferrer">this</a> and <a href="https://github.com/MatteoSchiavinato/Utilities/blob/master/SNP-dogg-v2.0" rel="nofollow noreferrer">this</a>. </p></li>
</ul>

CC BY-SA 4.0
2020-04-12T09:17:12.760Z
2020-04-12T09:17:12.760Z
schmat_90
<p>It looks like you mostly followed GATK?  It looks reasonable, though a lot will depend on how good your genome and gtf file are.</p>

CC BY-SA 4.0
2019-11-14T07:31:59.827Z
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.
Mohamed Moawad
CC BY-SA 4.0
2019-11-13T22:12:53.023Z
2019-11-13T22:12:53.023Z
swbarnes2
<p>Greetings,</p>

<p>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:</p>

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

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

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

<p>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).</p>

<p>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).</p>

<p>6- Calling variants again using BCFtools.</p>

<p>7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools. </p>

<p>8- Filtering variants using BCFtools.</p>

<p>9- Select variants SNPs and Indels in two separated files using GATK.</p>

<p>10- Predicting effects using SnpEff.</p>

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

<p>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? </p>

<p>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? </p>

<p>Your response is highly appreciated,</p>

<p>Kind regards,</p>

CC BY-SA 4.0
2019-11-13T13:44:25.617Z
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.
Matteo Ferla
CC BY-SA 4.0
2019-11-14T08:36:05.010Z
Do you mean blasting variants vcf file to the human genome? can you elaborate more?
Mohamed Moawad
CC BY-SA 4.0
2019-11-15T17:37:33.470Z
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](https://www.ensembl.org/info/data/ftp/index.html) 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.
Matteo Ferla
<snp><variant-calling><snpeff>
Questions regarding the detection of variants (SNPs and Indels) in mammal
# Questions regarding the detection of variants (SNPs and Indels) in mammal

(Asked by: username_0 on 2019-11-13)

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,

**Tags:** snp, variant-calling, snpeff
> **Comments:**
> **username_1**: 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.  
> **username_0**: Do you mean blasting variants vcf file to the human genome? can you elaborate more?  
> **username_1**: 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](https://www.ensembl.org/info/data/ftp/index.html) 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.  

---

## 2 answers

### Answer by username_2 on 2020-04-12. Score: 1

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](https://github.com/MatteoSchiavinato/Utilities/blob/master/add-read-groups.py).
* 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](https://github.com/MatteoSchiavinato/Utilities/blob/master/SNP-dogg) and [this](https://github.com/MatteoSchiavinato/Utilities/blob/master/SNP-dogg-v2.0).


---

### Answer by username_3 on 2019-11-13. Score: 0

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

> **Comments:**
> **username_0**: 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.  

---

