Sequence format: grouping sequences within a flat file

Question

Is there a sequence format allowing multiple sequences to be grouped across a flat file? The specific application is stacking separate nucleotide sequences against a given RNA secondary structure.

To rephrase the same question I'm looking to stack viral nucleotide data in different groups across a single flat file. FastaQ provides an example:

Traditional fastQ,

@sequence_name
AGCTGAT
+
...(quality score)

@sequence2
AGCTGAC
+
...(quality score)

What I'm looking for is a format,

@Gp1_virus1-2-3
AGACGTAG
+
AGAGACAA
+
AGAGACAT

@Gp2_virus1-2
AGATGTAG
+
AGATACAA

In the above example there are two groups comprising 3 (group 1) and 2 (group 2) individual sequences (viruses). Here 'group' represents different structures.

Stockholme format has potential. My understanding of this format, which is basic, is it stacks a sequence, together with structural information. What is sought is stacking multiple sequences and then a structure format. Any clarification would help.

There are multiple XML formats, which may be relevant.

Just to mention a 'flat file' is needed to allow manual inspection. Dumping JSON is a solution, but isn't friendly for optical inspection.

---

Application The immediate downstream application is structure/sequence viewers. I think at least one. viewer will accept Stockholme 1.0 format. rFAM is very similar (@terdon's comment)

Summarised Question The question could be rephrased as to whether multiple RNA families can be stacked inside one Stockholme 1.0 or else rFAM format?

Some of the members here are likely familiar with pFAM. So could multiple protein structures be stored and retrieved inside one pFAM file.

@gringer's contribution is very welcome.

Answer 1

An unmapped SAM file format allows for read groups. You can use samtools import to create unmapped SAM files from FASTQ files:

samtools import -R virus-1-2 in12.fq > 12_reads.sam
samtools import -R virus-1-2-3 in123.fq > 123_reads.sam
samtools merge out.sam 12_reads.sam 123_reads.sam

Answer 2

I'm awarding 'solved' to @gringer, because thats the context of the question.

The best solution I could find is stacking Stockholm (sth) formats together. In this example the RNA families are associated with one structure:

# STOCKHOLM 1.0
#=GF ID    UPSK
AF035635.1/619-641             UGAGUUCUCGAUCUCUAAAAUCG
M24804.1/82-104                UGAGUUCUCUAUCUCUAAAAUCG
J04373.1/6212-6234             UAAGUUCUCGAUCUUUAAAAUCG
M24803.1/1-23                  UAAGUUCUCGAUCUCUAAAAUCG
#=GC SS_cons                   .AAA....<<<<aaa....>>>>
//

Biopython allows sth format to be input.

from Bio import AlignIO
from pathlib import Path

file = 'test.sth'
path = '/Volumes/mypath/testSth'
alignment = AlignIO.read(open(Path(path, file)), "stockholm")
print("Alignment length %i" % alignment.get_alignment_length())
print (alignment.column_annotations['secondary_structure'])

Output,

Alignment length 23 
.AAA....<<<<aaa....>>>> # this is the RNA structure 

If we look the alignment object via

from pprint import pprint
pprint(vars(alignment)) # or via print(alignment.__dict__)

The output

{'_annotations': {}, '_per_col_annotations': {'secondary_structure': '.AAA....<<<<aaa....>>>>'},
'_records': [SeqRecord(seq=Seq('UGAGUUCUCGAUCUCUAAAAUCG'), id='AF035635.1/619-641', name='AF035635.1', description='AF035635.1/619-641', dbxrefs=[]), SeqRecord(seq=Seq('UGAGUUCUCUAUCUCUAAAAUCG'), id='M24804.1/82-104', name='M24804.1', description='M24804.1/82-104', dbxrefs=[]), SeqRecord(seq=Seq('UAAGUUCUCGAUCUUUAAAAUCG'), id='J04373.1/6212-6234', name='J04373.1', description='J04373.1/6212-6234', dbxrefs=[]), SeqRecord(seq=Seq('UAAGUUCUCGAUCUCUAAAAUCG'), id='M24803.1/1-23', name='M24803.1', description='M24803.1/1-23', dbxrefs=[])], 'annotations': {}} {'secondary_structure': '.AAA....<<<<aaa....>>>>'}

The reason vars(alignment) is required is to internally construct an alignment object in Stockholm format, here's the first part of the skeleton 'pseudo-code',

from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord

outAlign = MultipleSeqAlignment([])
record in alignment:
    name = record.id 
    myseq = record.seq
    sr = SeqRecord(Seq(name), myseq, '', '')
    outAlign.append(sr)

The above code constructs a basic alignment, the bit thats missing is to add the following key/value to the alignment object of {'secondary_structure': '.AAA....<<<<aaa....>>>>'}.

Its complicated, but you get use to it.

---

https://github.com/biopython/biopython/blob/master/Bio/AlignIO/StockholmIO.py