JBrowseMSA

An RNA family, from a model and six genomes

A riboswitch is a piece of an mRNA that folds up in front of the gene it controls and binds a small molecule through its fold. The sequence can drift as long as every base pair keeps two complementary sides, so a change on one side has to be met by a change on the other. This page takes one such family, finds its copies in six bacterial genomes, aligns them and opens the result with the base pairs drawn over the columns, where each row shows whether its pairs still match.

The pipeline runs curl, python3, cmsearch, cmalign and FastTree.

Prerequisites

Where the data comes from

Rfam release 15.1 for the family, NCBI RefSeq for the genomes and their annotation, and the RCSB PDB for one crystal structure.

The family

The SAM-I riboswitch, Rfam RF00162, sits in the leader of methionine and cysteine biosynthesis genes in Firmicutes and binds S-adenosylmethionine, the product those genes lead to. Rfam publishes two files for it: a seed alignment of 457 curated sequences carrying a consensus secondary structure, and a covariance model built from that seed. A covariance model scores a candidate on sequence and on fold at once, so it recognizes a copy that shares the fold and little of the sequence.

curl -o RF00162.cm https://rfam.org/family/RF00162/cm
curl -o RF00162.seed.sto https://rfam.org/family/RF00162/alignment/stockholm

1. Search six genomes

Six Firmicutes, downloaded one FASTA at a time from NCBI and concatenated: B. subtilis 168, Caldanaerobacter subterraneus subsp. tengcongensis MB4, Clostridium acetobutylicum ATCC 824, Listeria monocytogenes EGD-e, Staphylococcus aureus NCTC 8325 and Clostridioides difficile 630.

# --cut_ga applies the family's curated gathering threshold, the one Rfam uses
# to decide membership
cmsearch --cut_ga --noali --cpu 4 --tblout hits.tbl RF00162.cm all-genomes.fa

The search over six chromosomes takes seconds, and the table has 37 rows:

37 hits
  NC_000964.3  11 hits
  NC_003030.1  7 hits
  NC_003210.1  7 hits
  NC_003869.1  3 hits
  NC_007795.1  4 hits
  NC_009089.1  5 hits

B. subtilis has eleven, and it is the species in which Grundy and Henkin first described the S box regulon.

2. Name each row by the gene it leads

A riboswitch is a leader sequence, so each row is labeled with the gene downstream of it. NCBI serves the annotation of any window of a record as a feature table, so naming a hit takes one request. The hit at NC_000964.3:1180802-1180685 runs on the minus strand, so the window is the 1500 nt below it, read in the same direction:

curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore\
&id=NC_000964.3&seq_start=1179185&seq_stop=1180684&strand=2&rettype=ft&retmode=text"
90	>1500	gene
			gene	samT
			locus_tag	BSU_11010
90	>1500	CDS
			product	bifunctional homocysteine S-methyltransferase using (R,S)AdoMet and methylenetetrahydrofolate reductase [NAD(P)H]

The first CDS that reads forward in that frame is the gene downstream, and its symbol plus a four-letter species prefix is the row name: Bsub_samT. Where the annotation has no symbol the script uses the locus tag, as in Cdif_CD630_RS06510 in the figures below.

35 of 37 hits named by a gene, 2 with no same-strand gene within 1500 nt

Those two keep their coordinate instead, Cace_1073888 and Saur_2391667. Both sit in an intergenic gap whose flanking genes run the other way.

3. Align the hits to the model

cmalign aligns each hit back to the same model, so a column means one consensus position in every row:

cmalign --noprob -o cmalign.sto RF00162.cm hits.fa
37 rows x 187 columns, 108 consensus columns
62 columns are gaps in 90% or more of the rows

The model has 108 consensus positions and the alignment 187 columns. The other 79 are insertions, and most of them fall in one place: the hits run from 92 to 157 nt, and the rows at the long end put their extra sequence in the variable stem loop.

4. Put the pseudoknot back

A covariance model pairs columns the way brackets nest, and a pseudoknot is a pair that crosses a helix instead of nesting inside it. Rfam annotates the SAM-I pseudoknot in the seed’s SS_cons with letters, A for an opening end and a for its partner, and cmalign writes back only the nested helices its model holds.

Both alignments mark their consensus columns in #=GC RF, so the k-th consensus column of one is the k-th of the other. The build script copies the seed’s SS_cons across that correspondence, along with RNA_ligand_SAM, the columns Rfam marks as touching the ligand. The script leaves RNA_structural_elements behind, because the seed spells the element names across its own columns, insert columns included, and a per-column copy would keep the brackets and lose the labels.

4 pseudoknot pairs copied from the seed

5. Infer a tree

FastTree reads aligned DNA FASTA, so the build script converts the Stockholm with U written as T and Rfam’s insert gaps written as -:

FastTree -nt -gtr -nosupport tree-input.afa > family.nwk

The script then writes the Newick into the Stockholm as a #=GF NH line, where the viewer reads a tree from.

6. Open it

The alignment, its consensus structure, the SAM contacts and the tree now travel in one Stockholm file, so a ?data= link needs only its URL.

{
  "msaview": {
    "type": "MsaView",
    "colorSchemeName": "nucleotide",
    "colWidth": 8,
    "rowHeight": 15,
    "msaFilehandle": { "uri": "https://example.org/sam-riboswitch.sto" }
  }
}
The alignment holds 37 riboswitches over 187 columns. The tree on the left comes from the #=GF NH line, the Secondary-structure track from #=GC SS_cons, and the Base pairs arcs from that same string: blue for the nested helices, red for the pseudoknot. The other two tracks are the file’s remaining #=GC lines, drawn as they stand: RNA_ligand_SAM, with an X at each column that touches the ligand, and RF, the model’s own consensus sequence. The blank band across the middle is the variable stem loop, which holds most of the 62 columns that are gaps in 90% or more of the rows.

7. The pseudoknot

The same view with the pseudoknot marked. Its four pairs join columns 28-31 to 147-153, and the red arcs cross blue helix arcs instead of nesting with them. The structure track records the crossing with the letters A and a, and the arc track with the red color.

8. One helix, four base pairs

The view below zooms to the helix that runs from column 13 to column 46, where a change on one side of a pair comes with a matching change on the other.

Columns 18 and 41 in red, columns 35 and 36 in blue. Column 18 is U in 26 rows, C in 9 and G in 2; column 41 is G in 18, A in 17 and C in 2. Read as a pair they are UA in 17 rows, CG in 9, UG in 9 and GC in 2: four different base pairs, and a base pair in all 37 rows.

Columns 35 and 36 are the control, in the same frame. Column 36 varies more than column 18 does, its commonest base holding 57% of the rows against 26 of 37. SS_cons leaves both of them unpaired, and the last section searches the alignment for whatever does pair with them.

9. Where the ligand touches

Seven columns of the Rfam consensus are annotated as contacts with SAM, and the #=GC RNA_ligand_SAM track marks each of them with an X. Written into the snapshot as highlights as well, the same seven columns also get a labeled band down the rows.

Columns 7, 11, 50, 139, 140, 141 and 182, with the sequence logo track switched on above them. Averaged over the seven, the commonest base holds 99.2% of the rows. The paired columns average 75.5% and the columns that are neither paired nor in contact 78.0%.

Three of the seven, columns 139 to 141, are one side of a helix whose other side is at columns 48 to 50.

Columns 139, 140 and 141 read U, G and C down the rows. They pair with 50, 49 and 48, and those pairs are AU in 37 rows, CG in 37 and GC in 35. Unlike columns 18 and 41, each of these pairs holds one base pair in 35 to 37 of the 37 rows.

The aptamer whose crystal structure is PDB 2GIS, the structure Rfam took those contacts from, was cut out of one of the six genomes searched here. Comparing the construct against every row, the longest stretch it shares with one is 48 nt with Tten_TTE_RS08660, the riboswitch in front of the methylenetetrahydrofolate reductase gene of C. subterraneus, against 34 nt for the next row.

10. The tree

The tree with the whole alignment beside it at three pixels per column. Rows that lead the same gene in different species mostly do not come out together: the five rows named metK sit in three separate places. The tinted pair does come out together, and those two are the riboswitches of yxjG and yxjH, tandem paralogs whose leaders sit 1.4 kb apart on the B. subtilis chromosome. Four of the five C. difficile rows form another such clade.

The tree comes from 187 columns of a 108 nt RNA. It places related rows next to each other for reading, and its deeper branching is poorly supported.

11. Check it against the raw data

The figures show the consensus structure Rfam supplies. The last step of the build script re-derives what it can from the alignment alone, and prints three counts.

Every pair in SS_cons, over every row that has a base in both columns: 1223 complementary out of 1322, or 92.5%. Watson-Crick and GU both count, since GU is a pair a helix accepts.

paired columns: 1223/1322 = 92.5% can pair
pseudoknot only: 126/148 = 85.1% can pair
unpaired columns drawn at random: 8934/25415 = 35.2% can pair

The third line is the same statistic over pairs of columns drawn at random from the columns SS_cons leaves unpaired. Two nucleotide columns agree by chance often enough to matter, and 35.2% is how often.

The last check throws the structure away. For each of the ten most variable columns, search every other column for the one it pairs with best:

col  commonest base   best partner   SS_cons pairs them with
 58             31%     131 at  74%     131
133             32%      56 at  79%      56
131             34%      58 at  74%      58
 57             37%     132 at  91%     132
 56             40%     133 at  79%     133
132             40%      57 at  91%      57
  1             41%      11 at  73%     187
159             41%     172 at  86%     172
171             41%     146 at  70%     160
187             41%       5 at  81%       1
7 of 10 find the partner SS_cons names

Seven of the ten land on the column the published structure pairs them with, out of 98 candidates each. The three that miss are columns 1, 171 and 187, two of them at the ends of the alignment where fewer rows have a base at all.

The same search over the five most variable columns SS_cons leaves unpaired, which have no partner to find:

col  commonest base   best partner   SS_cons pairs them with
 60             44%     131 at  62%       -
136             44%      55 at  85%       -
165             49%      30 at  66%       -
 36             57%       5 at  76%       -
178             57%       5 at  84%       -

Column 36, the control from the helix figure, tops out at 76%. Two of the five reach the eighties, since a column scores well against any nearly invariant column carrying a complementary base. The seven hits above each found one specific column out of 98, the one SS_cons names, and these five did not.

Reproduce it end to end

curl -O https://raw.githubusercontent.com/GMOD/JBrowseMSA/main/docs/tutorials/scripts/build_rna_family.sh
bash build_rna_family.sh out/

The script fetches the model, the seed and the six genomes, runs every command above and writes out/sam-riboswitch.sto, the file the links on this page open. The numbers quoted here are the ones it prints.

See also

References