What is MGcount

MGcount is an RNA-seq quantification tool conceived to address ambiguous read alignments in a flexible way that is compatible with any biotype. This allows to extract more information from total RNA-seq datasets by simultaneous quantifying coding and non-coding transcripts, both small and long, where dealing with heterogeneous read-to-feature alignment ambiguities is key. When aligning reads to a reference genome, we distinguish between two types of ambiguous alignments:

• Multi-mappers: reads aligning to multiple genomic locations.
• Multi-overlaps: reads aligning to a genomic location with multiple annotated features.

To deal with the most frequent multi-overlaps, MGcount hierarchically assigns reads to small RNA, long RNA exons and long RNA introns, accounting for their length disparity. Subsequently, MGcount models read-to-feature alignments in a graph. This is exploited to detect and report expression of sequence-similar annotated loci, where reads systematically multi-map, as integrated features called communities.

Citation

Hita, A., Brocart, G., Fernandez, A., Rehmsmeier, M., Alemany, A. and Schvartzman, S., 2022. MGcount: a total RNA-seq quantification tool to address multi-mapping and multi-overlapping alignments ambiguity in non-coding transcripts. BMC Bioinformatics 23, 39 (2022). https://doi.org/10.1186/s12859-021-04544-3

System requirements

MGcount is built on the top of FeatureCounts, a well-known computational efficient software (Liao et. al., 2014). Please download it from the following link: http://subread.sourceforge.net/

Installation

MGcount is written in Python and is executed from the command line. You can either download the executable version (single binary file) or install it as a Python3 module.

pip3 install git+https://github.com/hitaandrea/MGcount.git

python3 -m mgcount [args]

chmod +x MGcount

MGcount [args]

1. Hierarchical assignation

To address the most frequent multi-overlapping ambiguous situations, reads are assigned to genomic annotated features in three pre-defined sequential rounds based on transcript body length: small RNA, long RNA exon and long RNA intron. MGcount prioritizes small RNA when a read aligns to a small RNA loci that is embedded within a long RNA. In the second round, alignments get assigned to long RNA exon features. In the final round, reads that haven’t been assigned to any small RNA or long RNA exon are assigned to introns if aligned within the full gene body length of a long RNA. Hence, all reads with at least one alignment overlapping with an annotated feature in the current round are assigned to such feature and skipped in subsequent rounds.

2. Multi-loci communities recognition

To quantify multi-mapping reads, MGcount builds a directed weighted graph {G=(V,E)} where each vertex (V) is an annotated feature and a pair of directional edges (E) connect two features for which common multi-mapping reads exist. Edge weights are defined as the ratio of multi-mapping reads between the two vertex normalized by the total number of reads assigned to the source vertex. Vertex weights are defined as the log-transformed number of assigned alignments (Fig c). Resultant graphs structures capture the multiple-loci topologies of different RNA biotypes. Two graphs are separately build for small-RNA and long-RNA features, using the full pool of input alignments. Subsequently, highly-related features are grouped together by minimizing the map equation with the communities detection approach described by Rosvall et al., 2008 .

3. Expression matrix generation

MGcount generates an output expression matrix for each hierarchical assignation round. These are appended together in a single output matrix. For each read, each alignment first gets an 1/N count, where N is the number of multi-mappers or residual multi-overlaps that survive the hierarchical assignment. Next, counts for annotated features which have been aggregated together in a community by the map equation are summed up. In this way, the systematic ambiguity in multi-mapping reads gets collapsed into a single MG community while the remaining signal is reported as fractional counts over distinct features.

Inputs

Three inputs are required to run the program:

• Input alignment file: a .txt file listing the paths to the .bam alignment input files by line.
• Annotations file: a .gtf file containing a set of RNA feature annotations
• Output directory path: a string specifying the path to the directory where MGcount outputs will be stored


( Please, use full-paths if you experience any problem)

Configurable arguments

Outputs

At the end of its execution, MGcount provides the following outputs:

• Count matrix: A matrix where each row corresponds to a feature as defined by feature_output (either single features or MG communities aggregating several features) and each column corresponds to one input BAM file.
• Features metadata: A table reporting: feature names matching row names in the count matrix, the counting round of hierarchical assignation, and its configuration parameters, a flag designing whether a feature belongs to an MG community, and the feature biotype.
• Multi-mapping graph adjacency matrix: A sparse adjacency matrix for each multi-mapping graph generated (small RNA and/or long RNA), stored as a symmetric, integer, squared matrix. Each matrix element stores the number of alignments that multi-map to a pair of features (defined by row and column), and the diagonal contains the total number of alignments per feature.
• Multi-mapping graph communities (MGcommunities): A table of MG communities linking each original feature in the GTF file with the resultant count matrix and metadata feature identifiers. It includes both unique features (which remain unmodified) and aggregated features (which are collapsed following MG communities). Also, the table stores the total number of alignments per feature.

T1 - Prepare inputs and execute

In the following tutorials, we use two sub-sampled human brain RNA-seq libraries as example to walk through MGcount execution. First of all, let’s create a folder and download the alignment .bam files of the two samples. (The downloading process might take a few minutes).

mkdir mgcount_tutorial
cd mgcount_tutorial

wget https://filedn.com/lTnUWxFTA93JTyX3Hvbdn2h/mgcount/tutorial_bamfiles.zip
unzip tutorial_bamfiles.zip -d input_bamfiles

To run MGcount, we need to provide the software with a .txt file specifying the paths to the input alignment files. Here, these are the two .bam alignment files we just downloaded. We can generate this file from the command line:

printf "$PWD/%s\n" input_bamfiles/* > input_bamfilenames.txt

The other required input is a .gtf file with transcript features annotations. MGcount repository provides with four ready .gtf files integrating annotations from several databases (see Hita et. al., 2022, Appendices, Methods, Database Integration). Next, we will download them and use the human annotations file for our execution example.

wget https://filedn.com/lTnUWxFTA93JTyX3Hvbdn2h/mgcount/integrated_annotations_gtf.zip
unzip integrated_annotations_gtf.zip -d annotations_gtf

Once we have both the .gtf and the input .bam files, we can simply run MGcount as an executable command-line program or as python3 module if installed via pip. A python module is run by calling python3 with the parameter -m and the name of the module, in this case “mgcount”. After this, we need to specify MGcount required arguments.

For this example, we parse the ready-integrated .gtf for the human genome, the .txt file we just created containing the path to the input alignment .bam files and a string designating the directory where MGcount outputs will be generated.

To reduce the computational time, we set the multicore parameter (-T) to “2” in order to parallelize the different steps of the algorithm by sample.

Run as an executable program:

MGcount -T 2 --gtf annotations_gtf/Homo_sapiens.GRCh38.gtf --outdir outputs --bam_infiles input_bamfilenames.txt

Run as a python3 module:

python3 -m mgcount -T 2 --gtf annotations_gtf/Homo_sapiens.GRCh38.gtf --outdir outputs --bam_infiles input_bamfilenames.txt

After MGcount run successfully finishes, your output directory should contain 6 new output files generate by MGcount including the RNA count matrix, the feature metadata table and two files containing the graph structure and the communities detected for long RNA and small RNA respectively.

Alternatively MGcount software might be invoked from an R console with the function “system”:

root_dir <- '~/mgcount_tutorial/'
system(paste0('MGcount -T 2',
              ' --gtf ',root_dir, 'annotations_gtf/Homo_sapiens.GRCh38.gtf',
              ' --outdir ' root_dir,'outputs ',
              ' --bam_infiles ',root_dir,'input_bamfilenames.txt'))

We are done!

T2 - Explore quantification outputs

In this tutorial, we will load the MGcount outputs we obtained in tutorial 1 and we will explore them from R. MGcount repository contains a few supporting scripts in R providing with side functionalities (manage annotations, visualize MGcount outputs…). It is possible to download the sole folder containing these scripts with the following shell command:

cd mgcount_tutorial
svn export https://github.com/hitaandrea/MGcount/trunk/R

Once the scripts are downloaded, we are ready to launch R. Please, start an R session and define the tutorial root directory as a variable.

root_dir <- '~/mgcount_tutorial/'

Further, this tutorial uses the following R packages. Please, install them with install.packages() if you wish to run this tutorial on your system.

library(dplyr)
library(Hmisc)
library(ggplot2)
library(ggpubr)
library(summarytools)

## Warning in fun(libname, pkgname): couldn't connect to display ":0"

source(paste0(root_dir,'R/integrate_gtf_annotations.R'))

The main output of MGcount is the count_matrix.csv file containing the feature by sample expression matrix. We can import it into R as any regular .csv.

counts <- read.csv(paste0(root_dir,'outputs/count_matrix.csv'), row.names = 1)
colnames(counts) <- sub('_Aligned.genome.dedup','',colnames(counts))

Each row in the count table is a feature for which expression has been quantified and each column is associated to a sample. By interrogating for the matrix dimension, we see the matrix contains two columns corresponding to the two human brain libraries. The matrix can be used as input for any RNA-seq downstream analysis.

dim(counts)

## [1] 46029     2

We can look at the total number of reads assigned from each library by summing up each of the two rows and compute the mean over the two libraries.

colSums(counts)

## Human_Brain_total_100ng_1_subsample Human_Brain_total_100ng_2_subsample 
##                             3260426                             3889868

print(paste('Mean counts:',mean(colSums(counts))))

## [1] "Mean counts: 3575147.255"

Let’s import now the feature_metadata output. This table reports feature-related attributes such as the assignation round to which the annotation belongs, a flag stating whether the feature is an aggregated community of annotations or an individual feature and the biotype associated to the feature. The feature identifiers in the counts_matrix and the feature_metadata match and therefore, this table can facilitate the extraction of particular features from the count matrix, e,g, a certain biotype, the exonic counts, the subset of features aggregated in communities, etc…)

feat_metadata <- read.csv(paste0(root_dir,'outputs/feature_metadata.csv'))

Here, for example, we profit from the feature_metadata table to look at counts distribution by assignation round.

df <- cbind('counts' = rowSums(counts), feat_metadata)
ggplot(df, aes(x = assignation_round, y = counts)) + geom_violin(fill = 'grey') + 
   scale_y_continuous(trans = 'log', lim = c(10,100000), breaks = c(10,100,1000,10000,100000)) + theme_pubclean() +
  geom_point(position = position_jitter(seed = 1, width = 0.4), size = 0.005, alpha = 0.2)

Below, we display a few random rows from the table for illustration.

feat_subset <- with(feat_metadata, feat_metadata[c(
   sample(which(assignation_round == 'small' & community_flag == "True"), 2),
   sample(which(assignation_round == 'small' & community_flag == "False"), 2),
   sample(which(assignation_round == 'long_exon' & community_flag == "True"), 1),
   sample(which(assignation_round == 'long_exon' & community_flag == "False"), 1),
   sample(which(assignation_round == 'long_intron' & community_flag == "True"), 1),
   sample(which(assignation_round == 'long_intron' & community_flag == "False"), 1)),])
kable(feat_subset) %>% scroll_box(height = "300px")

Next, we show how to generate a barplot showing the read distribution by biotype group. First lets load a table to group biotypes by category. We will use this to group less abundant biotypes in larger groups for visualization purposes.

bcats <- define_bcats()

feat_df <- merge(feat_metadata, bcats, by.x = 'feature_biotype', by.y = 'biotype', all.x = TRUE, all.y = FALSE) 

## Add exon/intron distinction to biogroup based on counting round
feat_df$biogroup <- as.character(feat_df$biogroup)
feat_df$biogroup[feat_df$biogroup == 'Protein_coding'] <- 'Protein_coding_exon'
feat_df$biogroup[feat_df$biogroup == 'Long_non_coding'] <- 'Long_non_coding_exon'  
feat_df$biogroup[feat_df$assignation_round == 'long_intron'] <-
  sub('exon','intron',feat_df$biogroup[feat_df$assignation_round == 'long_intron'])

feat_df$biogroup <- as.factor(feat_df$biogroup)
feat_df <- feat_df[match(rownames(counts), feat_df$feature),]

We then combine feature_metadata table with the count_matrix again and sum up the counts to get the total number of reads by biotype. We do this separately by biotype groups and small-non-coding individual biotypes to further represent the small non-coding spectrum.

## Generate biotype matrix
df <- data.frame(counts %>% group_by(feat_df$biogroup) %>% summarise_all(sum), 
                 check.names = FALSE); names(df)[1] <- 'biotype'
biotype <- reshape(df, idvar = 'biotype', varying = list(names(df)[-1]),
                      v.names = 'counts', times = names(df)[-1],
                      timevar = 'sn', direction = 'long')  
biotype$biotype <- as.character(biotype$biotype)
bgroups <- c("Hybrid","Long_pseudogenes","Protein_coding_intron","Long_non_coding_intron",
             "Protein_coding_exon","Long_non_coding_exon" ,"Short_non_coding","tRNA","rRNA")
biotype$biotype[biotype$biotype %nin% bgroups] <- 'Hybrid'
biotype$biotype <- factor(biotype$biotype, levels = bgroups)

## Generate non-coding biotype table
df <- data.frame(counts[feat_df$biocat == 'sNC',] %>%
                 group_by(feat_df$feature_biotype[feat_df$biocat == 'sNC']) %>% summarise_all(sum), 
                 check.names = FALSE); names(df)[1] <- 'biotype'
biotype_snc <- reshape(df, idvar = 'biotype', varying = list(names(df)[-1]),
               v.names = 'counts', times = names(df)[-1],
               timevar = 'sn', direction = 'long')  

Once we have extracted the counts by biotype group, let’s employ ggplot to visualize the read distribution profiles as barplots.

Reads distribution by biotype:

## ---- Abundance plot by biotype group
colP <- c('violetred1','slateblue1','darkgrey','lightgrey',
          'springgreen4','springgreen3','violetred4','tan3','tan4')
names(colP) <- bgroups
p1 <- ggplot(biotype, aes(x = sn, y = counts, group = biotype, fill = biotype)) +    
  geom_bar(stat = 'identity', colour = 'black', width = 0.8) + 
  coord_flip() + xlab('') + ylab('Number of assigned reads') + theme_pubclean() + 
  guides(fill=guide_legend(nrow=5)) + scale_fill_manual(values=colP) +  
  theme(legend.position = 'top', legend.title = element_blank())

## ---- Relative abundance plot by small non-coding biotype
colP <- c('#e6be97','#848CFF','#4c2382', '#50eb76','tomato','#FDFF87','#FFAE51','darkorchid1','gold1','sienna3')
p2 <- ggplot(biotype_snc, aes(x = sn, y = counts, group = biotype, fill = biotype)) +
  geom_bar(stat = 'identity', position = 'fill', colour = 'black', width = 0.8) + theme_pubclean() +
  coord_flip() + xlab('') + ylab('Proportion of small non-coding assigned reads') + 
  guides(fill=guide_legend(nrow=5)) +   scale_fill_manual(values = colP) +
  theme(axis.text.y = element_blank(), legend.position = 'top', legend.title = element_blank())

## Display plots together
ggarrange(p1, p2, ncol = 2, widths = c(2,1))

Finally, we import the multigraph_communities tables. These tables link each original feature in the .gtf with the resultant feature matching the count matrix and feature metadata identifiers. It includes both unique features (that remain unmodified) and aggregated features (that are collapsed following MG communities). Thus, we can track back the features grouped in each aggregated feature.

Here we will use the feature_metadata subset from before to explore the original features in the .gtf forming each new MG community feature.

csmall <- read.csv(paste0(root_dir,'outputs/multigraph_communities_small.csv'))
clong <- read.csv(paste0(root_dir,'outputs/multigraph_communities_long_exon.csv'))

kable(subset(csmall, feature %in% feat_subset$feature)) %>% scroll_box(height = "300px")

kable(subset(clong, feature %in% sub('-exon','',feat_subset$feature))) %>% scroll_box(height = "300px")

kable(subset(clong, feature %in% sub('-intron','',feat_subset$feature))) %>% scroll_box(height = "300px")

Also, we can exploit the multigraph_communities tables to extract stats from the MGcount run. Here we look at the proportion of aggregated features (community_flag variable) by biotype (feature_biotype). The output shows how small RNA biotypes tend to be more aggregated in communities because duplicated loci are more frequent.

kable(with(csmall, ctable(transcript_biotype, community_flag)))

kable(with(clong, ctable(gene_biotype, community_flag)))

We are done!

T3- Explore output multi-mapping graph

Below, we use a few functions defined in “mg_visualize.R” to graphically explore the multi-mapping graph topologies.

The function ‘mg_build’ takes the multi-mapping graph adjacency matrix and the MG communities and creates an igraph object that can be explored via igraph library. The next two functions (mg_plotset and mg_interactive) provide with the code to generate a few default plots given a Multi-Graph igraph object.

library(Matrix)
library(igraph)
library(plotly)
source(paste0(root_dir,'R/mg_visualize.R'))

As an example, we will load the small-RNA graph adjacency matrix and we will subset it to explore the sub-graph associated to microRNA features.

The adjacency matrix is stored as a sparse symmetric matrix in MatrixMarket format, which can be imported in R by the readMM function from ‘Matrix’ R package. We also import the table containing all annotated features and MG communities outputs.

dir.create(paste0(root_dir,'mgplots'))

## Import ml data and adjacency matrix for small assignation round
inputM <- readMM(paste0(root_dir,'outputs/multigraph_matrix_small.mtx'))
multiloci_table <- read.csv(paste0(root_dir,'outputs/multigraph_communities_small.csv'))

Next we subset the matrix by selecting the features under ‘miRNA category’. Each row in the communities table table corresponds to a feature annotation with non-zero assignments in the small round. The row index of each feature defines its column and row position in the multi-mapping graph adjacency matrix.

## Extract microRNA matrix subset
btype <- 'miRNA'
idx <- which(multiloci_table$transcript_biotype == btype)
inM <- inputM[idx,idx]; ml <- multiloci_table[idx,]

The necessary steps to convert the adjacency matrix and the communities table to an igraph object are in the mg_build function which outputs an igraph object.

## Generate microRNA graph object
attr <- 'transcript_name'
g = mg_build(inM, ml, attr)
g <- delete_vertices(g, V(g)[V(g)$weight < 1])

We can explore the patterns in the feature gene symbols established by the HGNC and link that to specific colors during visualization. Here we force annotations associated to a mature clipped 3-prime microRNA to be colored in orange and annotations associated to a mature clipped 5-prime microRNA to be colored in blue. For this, we look for “-3p” and “-5p” patterns in the transcript_name annotation symbol, which MGcount uses as default featue_output for small RNA.

## Extract microRNA mature extreme information from HUGO symol
V(g)$color_HUGO1 <- 'grey'
V(g)$color_HUGO1[grep('-3p',V(g)$feat)] <- 'orange'
V(g)$color_HUGO1[grep('-5p',V(g)$feat)] <- 'dodgerblue'

## Extract microRNA class from HUGO symbol
V(g)$color_HUGO2 <- 'grey'
idx <- grep('hsa-miR',as.character(V(g)$feat))
hugo <- as.factor(gsub('.*?([0-9]+).*', '\\1', V(g)$feat[idx]))
V(g)$color_HUGO2[idx] <- sample(grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)], 
                               length(unique(hugo)))[hugo]

The function mg_plotset generates a set of different plots of the multi-mapping graph and stores them as .png files with a user-given file path and prefix (plotfile argument).

We next use magick R package to load a few of the generated .png files into R. Each vertex is an annotated features with size proportional to its number of aligned reads. Each edge connects two features with shared multi-mappers with thickness proportional to the fraction of shared multi-mappers over the total alignments. Shared grey areas delineate MG communities. Vertices are colored according to the attribute “customColor”, which we just defined to be orange, blue and grey for “-3p”, “-5p” and absence of “-3p”/“-5p” patterns respectively. We may modify this as desired.

## Define plot colot
V(g)$customColor <- V(g)$color_HUGO1

## Standard plots
mg_plotset(g, plotfile= paste0(root_dir,'mgplots/mg_miRNA'))

Raw visualization of the graph

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA.png'))
plot(img)

Visualization of the graph colored by -3/-5 patterns in microRNA transcript symbol

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA_color.png'))
plot(img)

Visualization of the graph colored by -3/-5 patterns in microRNA transcript symbol and detected communities

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA_cl_color.png'))
plot(img)

Interactive visualization of the graph

To explore large graphs with interactive visualization tools such as zoom, we may use mg_interactive which creates a plotly object from the graph. Here the different communities are represented by colors.

mg_interactive(g, paste0(root_dir,'mgplots/mg_miRNA'))

We are done!