Alignment exercise answers - Revision history

Mick at 10:57, 19 December 2025

2025-12-19T10:57:43Z

← Older revision		Revision as of 12:57, 19 December 2025
Line 8:		Line 8:
	</ul>		</ul>

	<p><b>Q2.</b> ~~Four~~ header lines if you added the read group (<code>@RG</code>).</p>		<p><b>Q2.</b> Three header lines if you added the read group (<code>@RG</code>).</p>

	<p><b>Q3.</b> The genomic coordinate of the first read (<code>SRR8002634.1</code>) is: <b>position 166782</b>.</p>		<p><b>Q3.</b> The genomic coordinate of the first read (<code>SRR8002634.1</code>) is: <b>position 166782</b>.</p>
Line 19:		Line 19:

	<p><b>Q6.</b> Yes, <code>SRR8002634.11</code> is unmapped.		<p><b>Q6.</b> Yes, <code>SRR8002634.11</code> is unmapped.
	The SAM flag <code>0x4</code> (~~decimal~~ 4) indicates “read unmapped”.</p>		The SAM flag <code>0x4</code> (4) indicates “read unmapped”.</p>

	<p><b>Q7.</b> According to <code>flagstat</code>:</p>		<p><b>Q7.</b> According to <code>flagstat</code>:</p>
Line 27:		Line 27:
	</ul>		</ul>

	<p><b>Q8.</b> SAM was ~~467M~~ → BAM was ~~117M~~ → approximately <b>4× smaller</b>.</p>		<p><b>Q8.</b> SAM was 479M → BAM was 116M → approximately <b>4× smaller</b>.</p>

	<p><b>Q9.</b> BAM was ~~117M~~ → CRAM was ~~37M~~ → approximately <b>3.1× smaller</b>.</p>		<p><b>Q9.</b> BAM was 116M → CRAM was 38M → approximately <b>3.1× smaller</b>.</p>

	<p><b>Q10.</b></p>		<p><b>Q10.</b></p>
Line 87:		Line 87:

	<pre>samtools view NA19201.bam chr20 \| wc -l</pre>		<pre>samtools view NA19201.bam chr20 \| wc -l</pre>
	<p>→ ~3,~~449~~,~~006~~ (or very close, depending on how BWA breaks ties)</p>		<p>→ ~3,450,558 (or very close, depending on how BWA breaks ties)</p>

	<p>To exclude supplementary alignments (<code>0x800</code>):</p>		<p>To exclude supplementary alignments (<code>0x800</code>):</p>
	<pre>samtools view -F 0x800 NA19201.bam chr20 \| wc -l</pre>		<pre>samtools view -F 0x800 NA19201.bam chr20 \| wc -l</pre>
	<p>→ ~3,~~448~~,~~196~~</p>		<p>→ ~3,449,745</p>

	<p><b>Q19.</b> From the insert size distribution (via <code>plot-bamstats</code>):		<p><b>Q19.</b> From the insert size distribution (via <code>plot-bamstats</code>):
	Most common insert size is approximately <b>152 bp</b>.</p>		Most common insert size is approximately <b>152 bp</b>.</p>

	<p><b>Q20.</b> At the ~~IGV~~ position <code>chr20:35,581,362</code>:		<p><b>Q20.</b> At the position <code>chr20:35,581,362</code>:
	Reference base = <b>T</b>, sample base = <b>G</b>.</p>		Reference base = <b>T</b>, sample base = <b>G</b>.</p>

	<p><b>Q21.</b> Number of reads supporting the G allele: <b>8 reads</b>.</p>		<p><b>Q21.</b> Number of reads supporting the G allele: <b>8 reads</b> (the first line is the reference, the 2nd line is the sample consensus base). </p>

	<p><b>Q22.</b> ~~The variant lies in the gene: <b>FER1L4</b>~~.</p>		<p><b>Q22.</b> Homozygous as all reads support just a single allele.</p>

Mick at 12:23, 20 November 2025

2025-11-20T12:23:39Z

Show changes

WikiSysop: Created page with "'''Q1:''' 3 possible ways: * The file with the smaller file size contains the trimmed reads. * Peek in the file and determine which file contains reads of uneven lengths * Use fastqc, to determine which file contains overrepresented adapter sequences '''Q2:''' 4 lines if you have added the RG '''Q3:''' 166782 '''Q4:''' 0, this means that the probability of being mismapped is one. This means that this read cannot be confidently assigned to this position. '''Q5:''' The..."

2024-03-19T15:28:26Z

Created page with "'''Q1:''' 3 possible ways: * The file with the smaller file size contains the trimmed reads. * Peek in the file and determine which file contains reads of uneven lengths * Use fastqc, to determine which file contains overrepresented adapter sequences '''Q2:''' 4 lines if you have added the RG '''Q3:''' 166782 '''Q4:''' 0, this means that the probability of being mismapped is one. This means that this read cannot be confidently assigned to this position. '''Q5:''' The..."

New page

'''Q1:''' 3 possible ways:
* The file with the smaller file size contains the trimmed reads.
* Peek in the file and determine which file contains reads of uneven lengths
* Use fastqc, to determine which file contains overrepresented adapter sequences

'''Q2:''' 4 lines if you have added the RG

'''Q3:''' 166782

'''Q4:''' 0, this means that the probability of being mismapped is one. This means that this read cannot be confidently assigned to this position.

'''Q5:''' There are 3x16 and 5x0 so 3 on the - strand and 5 on the + strand

'''Q6:''' Yes it is unmapped, 4 means read unmapped.

'''Q7:''' 1183276 aligned out of 1252591 which is 94.47%. Therefore 100%-94.47% = 5.53%

'''Q8:''' Down from 467M to 117M so about 4X less.

'''Q9:''' Down from 117M to 37M so about 3.1X less.

'''Q10:''' Use:

<pre>
samtools view SRR8002634_1.sorted.bam NC_002516.2:2000000-3000000|wc -l
</pre>
155343

'''Q11:''' Use:

<pre>
samtools view -q 30 SRR8002634_1.sorted.bam NC_002516.2:2000000-3000000|wc -l
</pre>
149524

'''Q12:''' In SRR8002634_1.mosdepth.summary.txt we find:

<pre>
NC_002516.2 6264404 177098613 28.27 0 113
</pre>

The average is 28.27X, the maximum is 113 reads.

'''Q13:''' if you use the wrong genome reference, the number of aligned reads will likely go down. This is because as the phylogenetic distance increases, the more the genomes diverge, the less the reads are likely to land on the genome due to increased divergence. In the extreme, the reads aligned by serendipity. Conversely, if the reference genome used is closer phylogenetically to the species of the sample, the reads are more likely to align and the number of aligned reads will go up. For instance, 5689 (0.46%) aligned to the ''Yersinia pestis'' and 279217 (22.34%) aligned to the ''Pseudomonas aeruginosa''.

'''Q14''' the command line is:
<pre>
bwa mem -R "@RG\tID:RG26\tSM:YRB42" /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa /home/projects/22126_NGS/exercises/alignment/NA19201_1.fastq.gz
/home/projects/22126_NGS/exercises/alignment/NA19201_2.fastq.gz |samtools view -bS /dev/stdin | samtools sort /dev/stdin > NA19201.bam
</pre>

By the way, when we write:

<pre>
samtools view -bS /dev/stdin
</pre>

we can also write:

<pre>
samtools view -bS -
</pre>

It means the same thing and is quicker.

Very technically speaking, the command line above is not the most efficient. By default, BAM is a compressed binary format. Compressed here means it is run through '''gzip''' or more precisely '''bgzip''' (block version of '''gzip'''). It can also be uncompressed binary. The first call to '''samtools view''' transforms the data as compressed binary. Subsequently, '''samtools sort''' decompresses and re-compresses it. This is a bit wasteful, there is no point in compressing data if you're going to decompress it again in the subsequent step. Is smarter command line is:
<pre>
bwa mem -R "@RG\tID:RG26\tSM:YRB42" /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa /home/projects/22126_NGS/exercises/alignment/NA19201_1.fastq.gz
/home/projects/22126_NGS/exercises/alignment/NA19201_2.fastq.gz |samtools view -uS - | samtools sort /dev/stdin > NA19201.bam
</pre>

Note that we have changed '''samtools view -bS''' to '''samtools view -uS''', this says produced uncompressed BAM.

'''Q15''' 1) there are fewer intermediate files being produced 2) you guarantee parallel execution such that you are as fast as your bottleneck which in this case is mapping.

'''Q16''' Using:
<pre>
samtools flagstat NA19201.bam
</pre>
We find:
<pre>
3516557 + 0 mapped (99.83% : N/A)
</pre>
So 99.83%. Please note that this is unusually high as we knew that these are reads then mapped to chromosome 20.

'''Q17'''
<pre>
3490934 + 0 properly paired (99.23% : N/A)
</pre>

'''Q18''' Using:
<pre>
samtools view NA19201.bam chr20 |wc -l
3449006 (or 3448996)
</pre>
However, please note that you might get different numbers. Like mentioned in class, BWA tends to plate reads randomly in case of a tie. But expect the same ballpark. Also, to be more rigorous, as noted in the flagstat, there are 4620 supplementary alignments in total, so one would need to filter those out. Using

<pre>
samtools view -F0x800 NA19201.bam chr20 |wc -l
3448196 (or 3448186)
</pre>

'''Q19'''
<pre>
samtools stat NA19201.bam > NA19201.stat
plot-bamstats -p NA19201 NA19201.stat
</pre>

The most common insert size is 152.

'''Q20''':reference=T, sample=G

'''Q21''': 8 reads

'''Q22''': FER1L4