Jome0169
commented
3 years ago Anyways - with this all in mind, I'm interested in identifying the final set of mergers and splits, and seeing what proportion we agree/disagree on, and what we can say about these regions.
Pulling out Correct B73 Annotation Calls
I went ahead and dowloaded their file with the "Supported annotations according to our M2f procedure (Additonal File 10)" to my local computer directory here: /Users/feilab/Desktop/Monnahan_et_al/Monnahan_et_al_validated_merged_genes.txt
File looks like this:
The only way you're able to decipher the comparison being made here (this contains merge/split pairs from all combination of genomes) is by looking at the sequence ID name. B73 had the identifiers like: Zm00001d010484 where the 1d is the identification of B73.
With this in mind I'll need to filter all merged genes to those focusing on the B73 set. To do so, I ran a simple awk command looking for 1d in iether column 1 or column 2. Command" ❯ awk '$2 ~ /1d/ || $1 ~ /1d/ {print $0}' Monnahan_et_al_validated_merged_genes.txt | sort > Monnahan_et_al_validated_B73_genes.txt
When looking at the number of genes here you see something somewhat odd:
❯ awk '$2 ~ /1d/ || $1 ~ /1d/ {print $0}' Monnahan_et_al_validated_merged_genes.txt | sort | uniq | wc -l
14431443 isn' t what I was expecting. But this could be a trimmed down number to the 1383 mentioned in the manuscript. Quote:
Considering these split-genes along with the merged genes to which they corresponded, our analysis concerns 1275, 1383, and 2125 genes in the W22, B73, and PH207 annotations, respectively, corresponding to roughly 3–5% of all genes contained in these annotations.
Which would make sense because looking at the file - when pulling B73 specific annotations quite a few have the following values
Merged Splits M2f Call
Zm00004b039508 Zm00001d025384,Zm00001d025385 NA NA
Zm00004b039512 Zm00001d025395,Zm00001d025397 NA NA
Zm00004b040088 Zm00001d026086,Zm00001d026087 NA NA
Zm00004b040104 Zm00001d026105,Zm00001d026106 NA NA
Zm00004b040200 Zm00001d026220,Zm00001d026221 NA NATo get at the numbers of each annotation class they generated (NoCall, Merged, Split,NA) ran the following command and got the following values:
~/Desktop/Monnahan_et_al
❯ awk '{print $4}' Monnahan_et_al_validated_B73_genes.txt | sort | uniq -c
407 Merged
232 NA
552 NoCall
252 SplitInteresting... However this doesn't really tell us exactly how many of each merger/split class in B73 is actually supported. That "Merged" value above could correspond to a gene which is already Merged in B73 which is supported which is Split in Ph207. So I need to further process this.
Generate then the file of B73 true mergers as well as those splits in B73 which should stay as true splits, as well as those of No Calls and NAs Command:
awk '$2 ~ /1d/ && $4 == "Merged" {print $0}' Monnahan_et_al_validated_B73_genes.txt > Monnahan_et_al_validated_B73_genes.merged.txt
awk '$2 ~ /1d/ && $4 == "Split" {print $0}' Monnahan_et_al_validated_B73_genes.txt > Monnahan_et_al_validated_B73_genes.kept_split.txt
❯ awk '$2 ~ /1d/ && $4 == "NA" {print $0}' Monnahan_et_al_validated_B73_genes.txt > Monnahan_et_al_validated_B73_genes.NA.txt
❯ awk '$2 ~ /1d/ && $4 == "NoCall" {print $0}' Monnahan_et_al_validated_B73_genes.txt > Monnahan_et_al_validated_B73_genes.NoCall.txt
❯ wc -l Monnahan_et_al_validated_B73_genes.*
141 Monnahan_et_al_validated_B73_genes.NA.txt
240 Monnahan_et_al_validated_B73_genes.NoCall.txt
170 Monnahan_et_al_validated_B73_genes.kept_split.txt
96 Monnahan_et_al_validated_B73_genes.merged.txt
1443 Monnahan_et_al_validated_B73_genes.txt
2090 totalSo this is the set we actually want. This is the set of genes found in B73 which were thought to be needing either merger, or to remain a split. Excellent. From these we can start figuring out which of these intersect with my calls, and are maybe different, the same, or if my calls can potentially improve on their "No call" class.
Pulling My annotations
Generated my own set of passing genes using the following command: awk '{print $4}' merged_original.bed > Mendieta_et_al.passing.txt and put this in a sub directory named Mendieta. This is the list of all gene mergers which I found. I'll go ahead and grep these gene names later to see what proportion of identical genes is represented in this manuscript as compared to mine.
The manuscript Using multiple reference genomes to identify and resolve annotation inconsistencies is a recent publication by Monnahan et al from 2020. It details a clever approach which I appreciate which attempts to identify potentially split genes in the maize genome.
It does so by comparing multiple versions of maize genome annotations (B73 - Ph207 - W22) against one another, and compares them in a pairwise manner to identify genes with a "One gene hits many in another" genome approach. This can be somewhat seen in their figure 1.
The manuscript goes on to basically test whether these sets of genes should be merged - or kept as split genes. Asking the question - "which annotation is correct?" It does so by looking at RNA-seq data across multiple different tissue types which asking whether the 'split' genes involved are actually similarly expressed. This is summed in their metric M2f (Mean 2 fold expression difference) - Again figure 2 here:
It's definitely an interesting approach and one I appreciate. Something to note here though is that the cut off they're using - this "top 10 of m2f" metric is very stringent (which makes sense). Basically if the values of M2F for one of their calculated pairs doesn't fall on either side (the split or merged side) of the distribution, they do not "make a call" as it were. So this does mitigate conclusions they can make about some of the potential pairs they have. For instance - initially they identify 481 potential mergers in B73 - and end up identifying 96 potential mergers, 170 genes which should remain split, and 240 which were unable to be called (this adds up to 506??).