# How to compare sequences of genes obtained after whole-genome sequencing to reference genome?

My idea is to align whole-genome sequencing data (as fastq files after, 30× coverage, gDNA) to the reference mouse genome (NCBI), 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 (but it's closed), answer (gave a hint about whole reference mapping not only the locus of interest), question (not answered), question (useful I guess but not sure).

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 installed.

• 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? Oct 21 '21 at 11:53
• @terdon: Thanks for the comment. I've edited. But I didn't understand your question "What technique?". What do you mean?
– abc
Oct 21 '21 at 20:47
• 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 Bot
Oct 23 '21 at 14:40
• What have you tried? Nov 7 '21 at 20:10
• I've managed to get BAM and BAI files with usegalaxy.org. Now, I'm trying to visualize with IGV...
– abc
Nov 7 '21 at 22:33

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

# 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 - | \

# index the bamfile. makes a .bai file



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.

• conchoecia, thank you for your answer and suggestions. I did Bowtie2 at 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
Nov 7 '21 at 22:46
• bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT. Nov 8 '21 at 6:27
• @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 Nov 16 '21 at 5:49
• you can also try this command from the samtools docs "fastq to bam", 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 Nov 16 '21 at 5:56
• that all sounds expected :) Nov 17 '21 at 11:52

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