Question about m6A data processing and suppl table 13

[HenriettaHolze]

Hi,

Thank you for putting in so much effort into making this data available!
I have a question regarding the processing of the m6A methylation data.

I downloaded the supplementary table 13 to look at the methylation sites across transcripts.
I observed that for a genomic position, there is only ever 1 transcript associated.
It is hard to believe for me that if multiple transcripts are covering a genomic position, only 1 is ever methylated.
For some genes, there are multiple transcripts reported, and for some genomic position there are multiple genes. But never multiple transcripts of the same gene per genomic position.

Does this have something to do with how m6anet identifies m6a sites? Or is there some intermediate filtering step?
I'm trying to make sense of the script for fig. 6e-h but it uses python scripts and intermediate files that I cannot find.
https://github.com/GoekeLab/sg-nex-data/blob/master/manuscript/code/data%20analysis%20and%20visualization/Figure_6e-h.R

I would very much appreciate if you could clarify this.
If table 13 has some additional filtering, could you please make an unfiltered file available?

Lastly, could you please explain what the readcounts columns in the table mean exactly?
I observed different read counts for different genomic positions in a single transcript in a single sample. What does that mean? Does that indicate sequencing errors or incomplete reads?

Many thanks,
Henrietta

[HenriettaHolze]

I've now tried to reproduce table 13 by running the tutorial https://github.com/GoekeLab/sg-nex-data/blob/master/docs/SG-NEx_m6Anet_tutorial.md

I followed this tutorial to download the data.

And then installed m6anet with conda create --name m6anet -c bioconda m6anet -y
which installed m6anet 2.1.0

I updated the command to run the new version of m6anet:
m6anet inference --input_dir ./A549_directRNA_replicate6_run1 --out_dir ./A549_directRNA_replicate6_run1 --infer_mod_rate --n_processes 4

I went on to compare the results to the published data.

library(tidyverse)

Load suppl table 13.
Had to open and repair in excel, then save as csv.

iso_m6a <- read.csv("data/41592_2025_2623_MOESM16_ESM.csv")
dim(iso_m6a)

Load results from running m6anet.

data.site_proba <- read_csv("data/A549_directRNA_replicate6_run1/data.site_proba.csv")
head(data.site_proba)

Filter for sites with >= 50 reads, as was done in the publication.

data.site_proba <- data.site_proba %>%
  filter(n_reads >= 50)

Merge published data and m6anet results on transcript ID and position on transcript.

iso_m6a_data.site_proba_merged <-
  merge(
    iso_m6a,
    data.site_proba,
    by.x = c("tx_id", "transcriptomic_position"),
    by.y = c("transcript_id", "transcript_position"),
    all = T
  ) %>%
  select(
    gene_id,
    tx_id,
    transcriptomic_position,
    `A549_Rep6.Run1_readcounts`,
    `A549_Rep6.Run1_probability`,
    n_reads,
    probability_modified,
    mod_ratio
  ) %>%
  # remove rows that have no modification estimation for this sample
  filter(!is.na(A549_Rep6.Run1_probability) | !is.na(probability_modified))

dim(iso_m6a_data.site_proba_merged)

74820 rows.

sum(is.na(iso_m6a_data.site_proba_merged$A549_Rep6.Run1_probability))
sum(is.na(iso_m6a_data.site_proba_merged$probability_modified))

20577 positions have no data in suppl table 13 compared to m6anet results.
Only 2898 missing in m6anet output compared to published data.

iso_m6a_data.site_proba_merged %>%
  filter(probability_modified > 0.9) %>%
  filter(is.na(A549_Rep6.Run1_probability)) %>%
  nrow()

iso_m6a_data.site_proba_merged %>%
  filter(A549_Rep6.Run1_probability > 0.9) %>%
  filter(is.na(probability_modified)) %>%
  nrow()

No data in published table for 375 sites that are confidently modified.
Only 106 the other way round.

Compare the probability estimation.

iso_m6a_data.site_proba_merged %>%
  mutate(`A549_Rep6.Run1_probability` = replace_na(replace = -0.1, `A549_Rep6.Run1_probability`),
         probability_modified = replace_na(replace = -0.1, probability_modified)) %>%
  ggplot(aes(x = `A549_Rep6.Run1_probability`, y = probability_modified)) +
  geom_point(alpha = 0.1, size = 0.1) +
  theme_bw()

The probability of modification correlates but there are differences and as shown above many more sites in the re-processed data.

Compare read counts.

iso_m6a_data.site_proba_merged %>%
  mutate(`A549_Rep6.Run1_readcounts` = replace_na(replace = -0.1, `A549_Rep6.Run1_readcounts`),
         n_reads = replace_na(replace = -0.1, n_reads)) %>%
  ggplot(aes(x = A549_Rep6.Run1_readcounts, y = n_reads)) +
  geom_point(alpha = 0.1, size = 0.1) +
  theme_bw() +
  xlab("A549_Rep6.Run1_readcounts published") +
  ylab("n_reads m6anet re-run")

For 38662 sites the read count is exactly the same and for 12683 sites it differs (site reported in both tables).

Am I not using the same input data?
Has the published data been filtered?
Has the m6anet algorithm changed this much?

I also don't know how to get genomic positions for the novel sites that I get from re-running m6anet. Could you please make the script available? Thanks!

[yuukiiwa]

Hi Henrietta (@HenriettaHolze),

To answer your questions for your first comment above:

I downloaded the supplementary table 13 to look at the methylation sites across transcripts.
I observed that for a genomic position, there is only ever 1 transcript associated.
It is hard to believe for me that if multiple transcripts are covering a genomic position, only 1 is ever methylated.
For some genes, there are multiple transcripts reported, and for some genomic position there are multiple genes. But never multiple transcripts of the same gene per genomic position.
Table 13 is based on summarizing the probability modified at each specific genomic site on a chromosome, where multiple transcripts' probabilities modified at a site were averaged. We only included the transcript with the highest expression in the transcript_id column, which is why there is only one transcript reported for each position.

Does this have something to do with how m6anet identifies m6a sites? Or is there some intermediate filtering step?
I'm trying to make sense of the script for fig. 6e-h but it uses python scripts and intermediate files that I cannot find.
https://github.com/GoekeLab/sg-nex-data/blob/master/manuscript/code/data%20analysis%20and%20visualization/Figure_6e-h.R

I would very much appreciate if you could clarify this.
If table 13 has some additional filtering, could you please make an unfiltered file available?
Supplementary Table 13 is 0.9_avg_bycellline_matrix.csv in the script https://github.com/GoekeLab/sg-nex-data/blob/master/manuscript/code/data%20analysis%20and%20visualization/Figure_6e-h.R without filtering. The python scripts are for summarizing the bedgraph or m6anet inference outputs across the cell line replicates.

Lastly, could you please explain what the readcounts columns in the table mean exactly?
I observed different read counts for different genomic positions in a single transcript in a single sample. What does that mean? Does that indicate sequencing errors or incomplete reads?
The readcounts columns in the table show the number of reads supporting each specific genomic sites, which is summarized from the n_read columns across transcripts from the m6anet inference data.site_proba.csv output.

To answer the inconsistency seen in your second comment:

The analysis you did is based on one replicate from A549, which has two replicates, with transcript-level probability modified, whereas the analysis in the paper summairzed 1) all replicates within the same cell line and 2) all transcripts at the same genomic site. As you have mentioned, the inputs differ.

Functions to convert transcript to genomic position:

def readFasta(transcript_fasta,is_gff):
    fasta=open(transcript_fasta,"r")
    entries,separate_by_pipe="",False
    for ln in fasta:
        entries+=ln
    entries=entries.split(">")
    if len(entries[1].split("|"))>1:
        separate_by_pipe=True
    dict={}
    for entry in entries:
        entry=entry.split("\n")
        if len(entry[0].split()) > 0:
            id=entry[0].split('.')[0]
            seq="".join(entry[1:])
            dict[id]=[seq]
            if is_gff > 0:
                if separate_by_pipe == True:
                    g_id=info[1],split(".")[0]
                else:
                    g_id=entry[0].split("gene:")[1].split(".")[0]
                dict[id].append(g_id)
    return dict

def readAnnotation(gtf_or_gff):
    gtf=open(gtf_or_gff,"r")
    dict,is_gff={},0
    for ln in gtf:
        if not ln.startswith("#"):
            ln=ln.strip("\n").split("\t")
            if is_gff == 0:
                if ln[-1].startswith("ID") or ln[-1].startswith("Parent"):
                    is_gff = 1
                else:
                    is_gff = -1
            if is_gff < 0:
                if ln[2] == "transcript" or ln[2] == "exon":
                    chr,type,start,end=ln[0],ln[2],int(ln[3]),int(ln[4])
                    attrList=ln[-1].split(";")
                    attrDict={}
                    for k in attrList:
                        p=k.strip().split(" ")
                        if len(p) == 2:
                            attrDict[p[0]]=p[1].strip('\"')
                    ##tx_id=ln[-1].split('; transcript_id "')[1].split('";')[0]
                    ##g_id=ln[-1].split('gene_id "')[1].split('";')[0]
                    tx_id = attrDict["transcript_id"]
                    g_id = attrDict["gene_id"]
                    if tx_id not in dict:
                        dict[tx_id]={'chr':chr,'g_id':g_id,'strand':ln[6]}
                        if type not in dict[tx_id]:
                            if type == "transcript":
                                dict[tx_id][type]=(start,end)
                    else:
                        if type == "exon":
                            if type not in dict[tx_id]:
                                dict[tx_id][type]=[(start,end)]
                            else:
                                dict[tx_id][type].append((start,end))
            if is_gff > 0:
                if ln[2] == "exon" or ln[2] == "mRNA":
                    chr,type,start,end=ln[0],ln[2],int(ln[3]),int(ln[4])
                    tx_id=ln[-1].split('transcript:')[1].split(';')[0]
                    if ln[2] == "mRNA":
                        type="transcript"
                    if tx_id not in dict:
                        dict[tx_id]={'chr':chr,'strand':ln[6]}
                        if type == "transcript":
                            dict[tx_id][type]=(start,end)
                        if type == "exon":
                            dict[tx_id][type]=[(start,end)]
                    else:
                        if type == "transcript" and type not in dict[tx_id]:
                            dict[tx_id][type]=(start,end)
                        if type == "exon":
                            if type not in dict[tx_id]:
                                dict[tx_id][type]=[(start,end)]
                            else:
                                dict[tx_id][type].append((start,end))
    #convert genomic positions to tx positions
    if is_gff < 0:
        for id in dict:
            tx_pos,tx_start=[],0
            for pair in dict[id]["exon"]:
                tx_end=pair[1]-pair[0]+tx_start
                tx_pos.append((tx_start,tx_end))
                tx_start=tx_end+1
            dict[id]['tx_exon']=tx_pos
    else:
        for id in dict:
            tx_pos,tx_start=[],0
            if dict[id]["strand"] == "-":
                dict[id]["exon"].sort(key=lambda tup: tup[0], reverse=True)
            for pair in dict[id]["exon"]:
                tx_end=pair[1]-pair[0]+tx_start
                tx_pos.append((tx_start,tx_end))
                tx_start=tx_end+1
            dict[id]['tx_exon']=tx_pos
    return (dict,is_gff)

def t2g(gene_id,fasta_dict,annotation_dict,g2t_mapping,df_eventalign_index,readcount_min):
    tx_ids = []
    t2g_dict = {}
    transcripts = [tx for tx in annotation_dict if tx in g2t_mapping[gene_id]]
    n_reads = sum([len(df_eventalign_index.loc[tx]) for tx in transcripts])
    if n_reads >= readcount_min:
        for tx in transcripts:
            tx_seq = fasta_dict[tx][0]
            tx_contig = annotation_dict[tx]['chr']
            if tx_seq is None:
                continue
            for exon_num in range(len(annotation_dict[tx]['exon'])):
                g_interval=annotation_dict[tx]['exon'][exon_num]
                tx_interval=annotation_dict[tx]['tx_exon'][exon_num]
                for g_pos in range(g_interval[0],g_interval[1]+1): # Exclude the rims of exons.
                    dis_from_start = g_pos - g_interval[0]
                    if annotation_dict[tx]['strand'] == "+":
                        tx_pos = tx_interval[0] + dis_from_start
                    elif annotation_dict[tx]['strand'] == "-":
                        tx_pos = tx_interval[1] - dis_from_start
                    if (g_interval[0] <= g_pos < g_interval[0]+2) or (g_interval[1]-2 < g_pos <= g_interval[1]): # Todo: To improve the mapping
                        kmer = 'XXXXX'
                    else:
                        kmer = tx_seq[tx_pos-2:tx_pos+3]
                    t2g_dict[(tx,tx_pos)] = (tx_contig,gene_id,g_pos,kmer) # tx.contig is chromosome.
            tx_ids += [tx]

                
    return n_reads, tx_ids, t2g_dict

fasta used: https://zenodo.org/records/5707193/files/Homo_sapiens.GRCh38.cdna.ncrna_wtChrIs_modified.fa?download=1

gtf file used: https://ftp.ensembl.org/pub/release-91/gtf/homo_sapiens/Homo_sapiens.GRCh38.91.gtf.gz

Thanks!

Best wishes,
Yuk Kei