0
$\begingroup$

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

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

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

2 Answers 2

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
10
  • $\begingroup$ 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 $\endgroup$
    – abc
    Nov 7, 2021 at 22:46
  • 1
    $\begingroup$ bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT. $\endgroup$
    – conchoecia
    Nov 8, 2021 at 6:27
  • 1
    $\begingroup$ @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 $\endgroup$
    – Colin D
    Nov 16, 2021 at 5:49
  • 1
    $\begingroup$ 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 $\endgroup$
    – Colin D
    Nov 16, 2021 at 5:56
  • 1
    $\begingroup$ that all sounds expected :) $\endgroup$
    – Colin D
    Nov 17, 2021 at 11:52
1
$\begingroup$

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

Improve this answer
$\endgroup$
2

Your Answer

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

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