Summing position and read length exceeds chromosome size?

[IsabelMarleen]

Hello,

I have a question about the position and read lengths in the read headers produced by NanoSim. I am using the human_giab_hg002_sub1M_kitv14_dorado_v3.2.1 pretrained model to simulate sequences from hg002, hg003 and hg004. I observed some interesting behaviour when mapping the simulated reads and wanted to plot the amount of coverage the reads should provide based on where they are sampled from and compare that to what I get when I actually map the reads.

I thought I could use the position in the header as a start point, and the sum of the position and the length of the alignable middle region as the end point of where this read was obtained from in the original genome. When I try that, however, I get end positions for some of the reads that are larger than the chromosome they were sampled from. I can also see positions that are higher than the entire chromosome size.

Is this expected behaviour? If so, can you help me better understand what the position and read lengths mean?

I am also wondering about what the position means for reverse reads, is it the start or end point?

Below are the first few lines of some of my reads that show this pattern:

I created this in R using the code below. The input comb table was created by reading in all the headers from a nanosim file and processing it, and chr_sizes is another table with the lengths of the chromosomes in my reference genome (hg38).

comb %>%
left_join(chr_sizes) %>%
filter(pos+len_middle_region > size) %>%
mutate(len_diff = size-pos-len_middle_region) %>%
arrange(len_diff)

I'd really appreciate any insights you could give me. Thanks!

[lcoombe]

Hi @IsabelMarleen,

I just have a couple follow-up questions for you!

1. Could you provide one of the commands that you ran?
2. What version of NanoSim are you using?

Just want to make sure we are on the same page, and I'll look into this for you.

Thank you for your interest in NanoSim!
Lauren

[IsabelMarleen]

Hi @lcoombe ,

Thank you for your quick response.

The NanoSim command I ran to simulate the fragments is the following:
simulator.py genome -rg individual_genome.fa -c human_giab_hg002_sub1M_kitv14_dorado_v3.2.1 -n 2900000 --fastq -t 32 -min 401 --seed 4444600 -o sim_reads

My Nanosim version is bioconda::nanosim=3.2.2.

The individual genomes I am using to simulate from are created by running bcftools consensus -f $ref $vcf -o individual_genome.fa. The reference for this is hg38 (downloaded from https://ftp.ensembl.org/pub/release-114/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.{chromosome number}.fa.gz for each chromosome and then combined using cat) and the vcf files come from high-quality benchmarks (such as https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz). I use this individual genome to simulate the fragments that I then map against the original reference. Could the behaviour I am seeing be caused by a coordinate mismatch between the individual genome and the reference?

Thank you for having a look at this for me.

Best,
Isabel

[lcoombe]

Hi Isabel,

Thanks for the extra information!

I do think that the issue could be using hg38 instead of the individual_genome.fa file that you input to NanoSim. You would expect insertions/deletions to be applied during your bcftools step, which would skew the coordinates relative to hg38. As a sanity check, you could compile the same table as you did above, but using the chromosome lengths for your matched genome instead. If you ultimately still want the coordinates from hg38, you could try to compute offsets to liftover the coordinates from your individual genome to hg38, but I think that could lead to some issues (ie. what if the read started in an insertion region, etc). At that point, it honestly is probably easier to just re-map the reads to hg38.