alignment problem with reference genome

[Silvia-Giorgia-Signorini]

Hello Felix,
I’m a PhD student and I’m currently working on epigenetics data (DNA methylation) with limpets. I work in parallel with two different species, Patella caerulea and Patella ulyssiponensis, but I have some problems with the alignment of P. ulyssiponensis sequences to its reference genome. I know that it is a problem of the reference genome because I’ve aligned the sequence of P. ulyssiponensis to the other reference genome that I have (P. caerulea) and it worked well, even if the map efficiency was low (10.8%). The genome preparation was the same for the two species. And I also know that the script for the alignment is ok because it worked well with P. caerulea.

Final Alignment report

Sequence pairs analysed in total: 13359954
Number of paired-end alignments with a unique best hit: 5061646
Mapping efficiency: 37.9%
Sequence pairs with no alignments under any condition: 6812837
Sequence pairs did not map uniquely: 1485471
Sequence pairs which were discarded because genomic sequence could not be extracted: 5061646

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT: 0 ((converted) top strand)
GA/CT/CT: 0 (complementary to (converted) top strand)
GA/CT/GA: 0 (complementary to (converted) bottom strand)
CT/GA/GA: 2531542 ((converted) bottom strand)

Number of alignments to (merely theoretical) complementary strands being rejected in total: 0

Final Cytosine Methylation Report

Total number of C's analysed: 0

Total methylated C's in CpG context: 0
Total methylated C's in CHG context: 0
Total methylated C's in CHH context: 0
Total methylated C's in Unknown context: 0

Total unmethylated C's in CpG context: 0
Total unmethylated C's in CHG context: 0
Total unmethylated C's in CHH context: 0
Total unmethylated C's in Unknown context: 0

Can't determine percentage of methylated Cs in CpG context if value was 0
Can't determine percentage of methylated Cs in CHG context if value was 0
Can't determine percentage of methylated Cs in CHH context if value was 0
Can't determine percentage of methylated Cs in unknown context (CN or CHN) if value was 0

Bismark completed in 0d 3h 25m 51s

Here you can find other details:

• Reference genome: downloaded from NCBI (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_963678685.1/)
• Genome preparation script: bismark_genome_preparation --bowtie2 --verbose /cluster/scratch/ssignorini/testuly

Do you have any suggestion?
I really thank you in advance for your availability.
Silvia Signorini

[FelixKrueger]

Hi Silvia, something here definitely looks fishy (no pun intended). Could you provide me a with a few reads of your read library (ideally raw, untrimmed reads, and let me know which kind of library prep had been performed). Maybe 200K sequences (gzipped), as an email attachement?

I will in the meantime try to get the reference sequence and have it prepared. Cheers, Felix

[Silvia-Giorgia-Signorini]

Thanks a lot for your fast reply!
I have prepared individual libraries through the NEBNext Enzymatic Methyl-seq Kit (EM-seq™) and sequenced on an Illumina Novaseq 6000 sequencer and attached you can find a subsample of an individual (R1 and R2, raw data). Do you need more information?
Thanks a lot,
Silvia
subsampleR1_100K.fq.gz
subsampleR2_100K.fq.gz

[FelixKrueger]

I took your files, trimmed 10bp from both 5' ends, and ran alignments: (directional):

Sequence pairs analysed in total:       99872
Number of paired-end alignments with a unique best hit: 36615
Mapping efficiency:     36.7%

Sequence pairs with no alignments under any condition:  47409
Sequence pairs did not map uniquely:    15848
Sequence pairs which were discarded because genomic sequence could not be extracted:    0

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT:       18493   ((converted) top strand)
GA/CT/CT:       0       (complementary to (converted) top strand)
GA/CT/GA:       0       (complementary to (converted) bottom strand)
CT/GA/GA:       18122   ((converted) bottom strand)

Number of alignments to (merely theoretical) complementary strands being rejected in total:     0

Final Cytosine Methylation Report
=================================
Total number of C's analysed:   1769028

Total methylated C's in CpG context:    74903
Total methylated C's in CHG context:    1838
Total methylated C's in CHH context:    7128
Total methylated C's in Unknown context:        17

Total unmethylated C's in CpG context:  192829
Total unmethylated C's in CHG context:  296885
Total unmethylated C's in CHH context:  1195445
Total unmethylated C's in Unknown context:      430

C methylated in CpG context:    28.0%
C methylated in CHG context:    0.6%
C methylated in CHH context:    0.6%
C methylated in unknown context (CN or CHN):    3.8%


Bismark completed in 0d 0h 1m 47s

So all works well, there is a 50/50 split between top and bottom strand. I am not exactly sure what went wrong on your side, but something did... Do you have enough resources available at all?

[Silvia-Giorgia-Signorini]

I apologize for my late reply but I was trying to figure out what was my problem. I think I had an issue with the original fasta file of the reference genome, now I fixed it thanks to your indications.
I really thank you for your time!
Best wishes
Silvia

[FelixKrueger]

That's good to hear! All the best