<p>Surprisingly there aren't actually many helpful one-liners to map Illumina WGS reads on the internet - everyone just sort of assumes that people can figure it out from the manual. Many people end up doing it in a not-so-ideal way, in that they stop at the <code>.sam</code> files that <code>bwa mem</code> makes that take up a lot of disk space.</p>
<p>I like this incantation:</p>
<pre><code># index the reference
bwa index reference_genome.fasta

# map your reads to the reference (bwa mem)
# compress that output to a bam file (samtools view -hb)
# sort the bam file (samtools sort)

bwa mem -t &lt;number_of_threads&gt; reference_genome.fasta your_sample.R1.fastq.gz your_sample.R2.fastq.gz | \
samtools view -hb - | \
samtools sort - &gt; your_reads_to_reference.sorted.bam

# index the bamfile. makes a .bai file

samtools index your_reads_to_reference.sorted.bam
</code></pre>
<p>With this sorted <code>.bam</code>, the <code>.bai</code>, and the reference genome you can now use <a href="https://software.broadinstitute.org/software/igv/" rel="nofollow noreferrer">IGV</a> to directly look at the reads mapped to your genome.</p>
<p>Also with the bam, you can generate a VCF file using <a href="https://github.com/freebayes/freebayes" rel="nofollow noreferrer"><code>freebayes</code></a> or another tool like <a href="https://www.biostars.org/p/63429/" rel="nofollow noreferrer"><code>samtools mpileup</code></a>. This will give you a file of how your sample is different from the reference genome. You can use <a href="https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html" rel="nofollow noreferrer"><code>bedtools intersect</code></a> to get the variants only from your gene/regions of choice.</p>

CC BY-SA 4.0
2021-11-07T22:46:28.447Z
conchoecia, thank you for your answer and suggestions. I did Bowtie2 at https://usegalaxy.org. Liam McIntyre and you suggest `bwa mem`, `minimap2`, or `ngmlr`. My reads are ~150nt. What tool should I choose? I guess for 150nt `bwa mem` is better. Currently with `.bam` and `.bai` files I'm trying to visualize data in IGV, but yes, genes of interest are preferred to show. I hope to completely switch to command line
abc
CC BY-SA 4.0
2021-11-08T06:27:42.390Z
bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT.
conchoecia
CC BY-SA 4.0
2021-11-08T06:53:13.630Z
Then I'll do bwa mem. Thank you
abc
CC BY-SA 4.0
2021-11-11T22:27:08.973Z
according to timpo0's comment https://github.com/jts/methylation-analysis/issues/13#issuecomment-370628831 with reference to  https://github.com/lh3/bwa ("For all the algorithms, BWA first needs to construct the FM-index for the reference genome (the index command)."): `bwa index reference.fasta` precede the rest lines
abc
CC BY-SA 4.0
2021-11-13T07:21:19.890Z
.bam file's size is 0 bytes...
abc
CC BY-SA 4.0
2021-11-16T05:49:36.757Z
@abc I also got 0 bytes also.  here is a possible helper `bwa index ref.fa; bwa mem ref.fa out1.fq out2.fq |  samtools sort - -o out.sorted.bam`. I am using samtools 1.14 and bwa-0.7.17-r1188
Colin D
CC BY-SA 4.0
2021-11-16T05:56:15.430Z
you can also try this command from the samtools docs "fastq to bam", http://www.htslib.org/workflow/#fastq_to_bam just swap `minimap2` with `bwa mem` 
 e.g. `bwa mem ref.fa out1.fq out2.fq | samtools fixmate -u -m - - | samtools sort -u -@2 - | samtools markdup -@8 --reference ref.fa - final.cram`
Colin D
CC BY-SA 4.0
2021-11-17T07:09:06.700Z
@ColinD: thanks for suggestion, 9h run (1st code) resulted in 12 out.sorted.bam.tmp.00XX.bam files of ~260Mb, the run is still going. A file is generated ~ every 30mins. In shell I see `[M::process] read 70190 sequences (10000248 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (8, 26835, 0, 4)
[] skip orientation FF as there are not enough pairs
[] analyzing insert size distribution for orientation FR...
[] (25, 50, 75) percentile: (172, 245, 331)
[] low and high boundaries for computing mean and std.dev: (1, 649)
[] mean and std.dev: (257.71, 118.14) <...>`
abc
CC BY-SA 4.0
2021-11-17T07:11:48.210Z
`<...> [] low and high boundaries for proper pairs: (1, 808)
[] skip orientation RF as there are not enough pairs
[] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 70060 reads in 64.900 CPU sec, 64.655 real sec`
abc
CC BY-SA 4.0
2021-11-17T11:52:17.473Z
that all sounds expected :)
Colin D
CC BY-SA 4.0
2021-11-07T17:03:58.743Z
2021-11-21T00:50:50.803Z
2021-11-21T00:50:50.803Z
conchoecia
<p>This is very easy. First map all reads to the reference sequence for the organism. If they are short reads use BWA mem. If long reads use minimap2 or ngmlr. Once mapped view the locus of interest in igv (GUI). If you need to extract the reads at the locus pass the coordinates to samtools view (cli).</p>

CC BY-SA 4.0
2021-10-22T06:19:22.917Z
MvIntyre: thank you for the answer, may be you can give advice any source on the topic as well ? I found BWA-MEM: https://github.com/lh3/bwa and https://www.ridom.de/u/BWA_Mapper.html, IGV: http://software.broadinstitute.org/software/igv/home, Samtools: http://www.htslib.org/
abc
CC BY-SA 4.0
2021-10-23T10:13:58.123Z
They are all git projects
Liam McIntyre
CC BY-SA 4.0
2021-10-21T20:39:10.097Z
2021-10-21T20:39:10.097Z
Liam McIntyre
<p>My idea is to align whole-genome sequencing data (as fastq files after, 30× coverage, gDNA) to the reference mouse genome (<a href="https://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.27" rel="nofollow noreferrer">NCBI</a>), extract the immunoglobulin locus and compare it to the reference. I think the question is too broad but maybe a hint how to start ? (I'm new to the field.)</p>
<p>Task is: to check for difference between reference genome and a sample for 100-200 genes.</p>
<p>I checked these: <a href="https://bioinformatics.stackexchange.com/q/16303/14406">question</a> (but it's closed), <a href="https://bioinformatics.stackexchange.com/a/11396/14406">answer</a> (gave a hint about whole reference mapping not only the locus of interest), <a href="https://bioinformatics.stackexchange.com/q/10633/14406">question</a> (not answered), <a href="https://bioinformatics.stackexchange.com/q/10633/14406">question</a> (useful I guess but not sure).</p>
<p>Additional info:
Instrument Platform: ILLUMINA<br />
Instrument Model: HiSeq X Ten<br />
Library Layout: PAIRED<br />
Library Strategy: WGS<br />
Library Source: GENOMIC</p>
<p>Typical sequence length distribution from FastQC: 150nt.</p>
<p>I have Win10 with <a href="https://bioinformatics.stackexchange.com/q/13245">WSL 1</a> installed.</p>

CC BY-SA 4.0
2021-10-21T11:53:01.987Z
Please [edit] your question and give us more details. What species are you working with? What technique? Is this DNAseq or RNAseq or something else? Do you know that the locus is or is that what you want to find out?
terdon
CC BY-SA 4.0
2021-10-21T20:47:55.367Z
@terdon: Thanks for the comment. I've edited. But I didn't understand your question "What technique?". What do you mean?
abc
CC BY-SA 4.0
2021-10-23T14:40:17.490Z
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Community
CC BY-SA 4.0
2021-11-07T20:10:18.057Z
What have you tried?
swbarnes2
CC BY-SA 4.0
2021-11-07T22:33:16.437Z
I've managed to get BAM and BAI files with https://usegalaxy.org. Now, I'm trying to visualize with IGV...
abc
<phylogenetics><fastq><reference-genome>
How to compare sequences of genes obtained after whole-genome sequencing to reference genome?
# How to compare sequences of genes obtained after whole-genome sequencing to reference genome?

(Asked by: username_0 on 2021-10-21)

My idea is to align whole-genome sequencing data (as fastq files after, 30× coverage, gDNA) to the reference mouse genome ([NCBI](https://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.27)), extract the immunoglobulin locus and compare it to the reference. I think the question is too broad but maybe a hint how to start ? (I'm new to the field.)

Task is: to check for difference between reference genome and a sample for 100-200 genes.

I checked these: [question](https://bioinformatics.stackexchange.com/q/16303/14406) (but it's closed), [answer](https://bioinformatics.stackexchange.com/a/11396/14406) (gave a hint about whole reference mapping not only the locus of interest), [question](https://bioinformatics.stackexchange.com/q/10633/14406) (not answered), [question](https://bioinformatics.stackexchange.com/q/10633/14406) (useful I guess but not sure).

Additional info:
Instrument Platform: ILLUMINA  
Instrument Model: HiSeq X Ten  
Library Layout: PAIRED  
Library Strategy: WGS  
Library Source: GENOMIC

Typical sequence length distribution from FastQC: 150nt.

I have Win10 with [WSL 1](https://bioinformatics.stackexchange.com/q/13245) installed.

**Tags:** phylogenetics, fastq, reference-genome
> **Comments:**
> **username_1**: Please [edit] your question and give us more details. What species are you working with? What technique? Is this DNAseq or RNAseq or something else? Do you know that the locus is or is that what you want to find out?  
> **username_0**: @terdon: Thanks for the comment. I've edited. But I didn't understand your question "What technique?". What do you mean?  
> **username_2**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
> **username_3**: What have you tried?  
> **username_0**: I've managed to get BAM and BAI files with https://usegalaxy.org. Now, I'm trying to visualize with IGV...  

---

## 2 answers

### Answer by username_4 on 2021-11-07. Score: 2

Surprisingly there aren't actually many helpful one-liners to map Illumina WGS reads on the internet - everyone just sort of assumes that people can figure it out from the manual. Many people end up doing it in a not-so-ideal way, in that they stop at the .sam files that bwa mem makes that take up a lot of disk space.
I like this incantation:
# index the reference
bwa index reference_genome.fasta

# map your reads to the reference (bwa mem)
# compress that output to a bam file (samtools view -hb)
# sort the bam file (samtools sort)

bwa mem -t  reference_genome.fasta your_sample.R1.fastq.gz your_sample.R2.fastq.gz | \
samtools view -hb - | \
samtools sort - > your_reads_to_reference.sorted.bam

# index the bamfile. makes a .bai file

samtools index your_reads_to_reference.sorted.bam

With this sorted .bam, the .bai, and the reference genome you can now use IGV to directly look at the reads mapped to your genome.
Also with the bam, you can generate a VCF file using freebayes or another tool like samtools mpileup. This will give you a file of how your sample is different from the reference genome. You can use bedtools intersect to get the variants only from your gene/regions of choice.

> **Comments:**
> **username_0**: conchoecia, thank you for your answer and suggestions. I did Bowtie2 at https://usegalaxy.org. Liam McIntyre and you suggest `bwa mem`, `minimap2`, or `ngmlr`. My reads are ~150nt. What tool should I choose? I guess for 150nt `bwa mem` is better. Currently with `.bam` and `.bai` files I'm trying to visualize data in IGV, but yes, genes of interest are preferred to show. I hope to completely switch to command line  
> **username_4**: bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT.  
> **username_0**: Then I'll do bwa mem. Thank you  
> **username_0**: according to timpo0's comment https://github.com/jts/methylation-analysis/issues/13#issuecomment-370628831 with reference to  https://github.com/lh3/bwa ("For all the algorithms, BWA first needs to construct the FM-index for the reference genome (the index command)."): `bwa index reference.fasta` precede the rest lines  
> **username_0**: .bam file's size is 0 bytes...  
> **username_5**: @abc I also got 0 bytes also.  here is a possible helper `bwa index ref.fa; bwa mem ref.fa out1.fq out2.fq |  samtools sort - -o out.sorted.bam`. I am using samtools 1.14 and bwa-0.7.17-r1188  
> **username_5**: you can also try this command from the samtools docs "fastq to bam", http://www.htslib.org/workflow/#fastq_to_bam just swap `minimap2` with `bwa mem`   e.g. `bwa mem ref.fa out1.fq out2.fq | samtools fixmate -u -m - - | samtools sort -u -@2 - | samtools markdup -@8 --reference ref.fa - final.cram`  
> **username_0**: @ColinD: thanks for suggestion, 9h run (1st code) resulted in 12 out.sorted.bam.tmp.00XX.bam files of ~260Mb, the run is still going. A file is generated ~ every 30mins. In shell I see `[M::process] read 70190 sequences (10000248 bp)... [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (8, 26835, 0, 4) [] skip orientation FF as there are not enough pairs [] analyzing insert size distribution for orientation FR... [] (25, 50, 75) percentile: (172, 245, 331) [] low and high boundaries for computing mean and std.dev: (1, 649) [] mean and std.dev: (257.71, 118.14) <...>`  
> **username_0**: `<...> [] low and high boundaries for proper pairs: (1, 808) [] skip orientation RF as there are not enough pairs [] skip orientation RR as there are not enough pairs [M::mem_process_seqs] Processed 70060 reads in 64.900 CPU sec, 64.655 real sec`  
> **username_5**: that all sounds expected :)  

---

### Answer by username_6 on 2021-10-21. Score: 1

This is very easy. First map all reads to the reference sequence for the organism. If they are short reads use BWA mem. If long reads use minimap2 or ngmlr. Once mapped view the locus of interest in igv (GUI). If you need to extract the reads at the locus pass the coordinates to samtools view (cli).

> **Comments:**
> **username_0**: MvIntyre: thank you for the answer, may be you can give advice any source on the topic as well ? I found BWA-MEM: https://github.com/lh3/bwa and https://www.ridom.de/u/BWA_Mapper.html, IGV: http://software.broadinstitute.org/software/igv/home, Samtools: http://www.htslib.org/  
> **username_6**: They are all git projects  

---

