Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-07T22:46:28.447Z
Answers.Text: 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
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-08T06:27:42.390Z
Answers.Text: bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-08T06:53:13.630Z
Answers.Text: Then I'll do bwa mem. Thank you
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-11T22:27:08.973Z
Answers.Text: 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
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-13T07:21:19.890Z
Answers.Text: .bam file's size is 0 bytes...
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-16T05:49:36.757Z
Answers.Text: @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
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-16T05:56:15.430Z
Answers.Text: 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`
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-17T07:09:06.700Z
Answers.Text: @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) <...>`
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-17T07:11:48.210Z
Answers.Text: `<...> [] 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`
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-17T11:52:17.473Z
Answers.Text: that all sounds expected :)
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-07T17:03:58.743Z
Answers.LastActivityDate: 2021-11-21T00:50:50.803Z
Answers.LastEditDate: 2021-11-21T00:50:50.803Z
Answers.UserDisplayName: conchoecia
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-10-22T06:19:22.917Z
Answers.Text: 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/
Answers.UserDisplayName: abc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-10-23T10:13:58.123Z
Answers.Text: They are all git projects
Answers.UserDisplayName: Liam McIntyre
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-10-21T20:39:10.097Z
Answers.LastActivityDate: 2021-10-21T20:39:10.097Z
Answers.UserDisplayName: Liam McIntyre
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-10-21T11:53:01.987Z
Comments.Text: 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?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-10-21T20:47:55.367Z
Comments.Text: @terdon: Thanks for the comment. I've edited. But I didn't understand your question "What technique?". What do you mean?
Comments.UserDisplayName: abc
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-10-23T14:40:17.490Z
Comments.Text: 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.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-07T20:10:18.057Z
Comments.Text: What have you tried?
Comments.UserDisplayName: swbarnes2
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-07T22:33:16.437Z
Comments.Text: I've managed to get BAM and BAI files with https://usegalaxy.org. Now, I'm trying to visualize with IGV...
Comments.UserDisplayName: abc
Tags: <phylogenetics><fastq><reference-genome>
Title: How to compare sequences of genes obtained after whole-genome sequencing to reference genome?
ThreadText: # 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**: @username_1: 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 <number_of_threads> 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](https://software.broadinstitute.org/software/igv/) to directly look at the reads mapped to your genome.
  
  Also with the bam, you can generate a VCF file using [`freebayes`](https://github.com/freebayes/freebayes) or another tool like [`samtools mpileup`](https://www.biostars.org/p/63429/). This will give you a file of how your sample is different from the reference genome. You can use [`bedtools intersect`](https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html) to get the variants only from your gene/regions of choice.
  
  > **Comments:**
  > **username_0**: username_4, thank you for your answer and suggestions. I did Bowtie2 at https://usegalaxy.org. username_6 username_6 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**: @username_0 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**: @username_5: 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  
  
  ---
  
  