Hello everyone,
I have a multiple query genomes aligned to a target genome stored in a CRAM file where each query genome is represented by a read group. I want to extract a multiple sequence alignment for a specific region, say for instance chr1:1,000,000-1,111,000. Usual output formats such as MAF, ALN or aligned FASTA would be fine.
I tried samtools view but it reports all sequences which overlap the window, which can be extremely long in case there are two highly related genomes in the dataset, and I did not find a way or a tool to clip the output to the region only. Also, I do not know how to convert from SAM/BAM/CRAM to aligned FASTA (that is, FASTA with gaps) or other multiple alignment formats.
I also tried samtools tview, and its option -p allows me to jump directly to the region of interest in interactive mode. Unfortunately I have not found a way to export its output in any useful format. In addition, there seems to be a bug that prevented the --reference option to provide the target genome sequence.
I know that storing genome alignments in CRAM format is unconventional, but this saves a lot of space and the ability index the multiple genome alignment and extract regions appears powerful… if the output format would be useful.
Do you know a solution?
Charles
•
written
by
Charles Plessy
1 answer
I wrote sam2logo https://jvarkit.readthedocs.io/en/latest/SAM4WebLogo/
for 111,000 bp ? it should work
Thanks Pierre and GenoMax for your suggestions. I could not manage to run sam2fasta but was able to run SAM4WebLogo and gofasta. Unfortunately, I had problems in both cases: SAM4WebLogo did not cut the sequences at the right position, and gofasta missed some sequences and did not output inserted bases. I attached a screenshot. The data is from an alignment that does not use public data, but I can share it in private if there is interest. Here are my command lines:
samtools tview --reference ref.bgzip.fa.gz aln.cram XSR:1111431-1111661 -p XSR:1111431
./jvarkit.jar sam4weblogo --reference ref.bgzip.fa.gz -r XSR:1111431-1111641 aln.cram > sam4weblogo.fasta
samtools view -h --reference ref.bgzip.fa.gz aln.cram XSR:1111431-1111641 | ./gofasta-linux-amd64 sam toMultiAlign -r ref.bgzip.fa.gz -w 80 --start 1111431 --end 1111641 > gofasta.fasta
Log in to answer this question.
-
What is wrong with my NCBI esearch command for 16S rRNA?
written by Morgan S. •
-
Comparing Alignment Files (CRAM)
written by joe_genome •
-
How to convert exonerate C4 alignment output to gff?
written by Junior •
-
Tool: A pairwise genome alignment pipeline using LAST and Nextflow
written by Charles Plessy
-
Software for Unaligning Sequences
written by printingdots •
-
What is the most suitable file format for storing sequence alignments and associated annotations?
written by Joseph Hughes
-
Can I get Multiple Sequence Alignments from Satsuma Alignments
written by ben.ward •
-
covert STOCKHOLM alignment to FASTA alignment
written by venu
-
Reference independent bacterial genome visualization
written by h.botond •
-
Psiblast Fasta Formatted Output
written by Reyhaneh
It sounds like you're trying to re-implement pangenomics using linear aligners.
Have you considered using minigraph-cactus and vg/odgi downstream instead ? I think you'll get more joy out of this for most downstream applications.
Thanks for the suggestion but I am working at large evolutionary distances (for instance aligning the human genome to 75 primate chromosome assemblies ranging from lemurs to great apes) and if I understand well, pangenome software are not designed for that. Also the pipeline I use (https://nf-co.re/pairgenomealign/) completes all the alignments in parallel within a hour, and to the best of my knowledge, cactus has the reputation of being very slow…
progressiveCactus is built for large evolutionary distances https://www.nature.com/articles/s41586-020-2871-y (https://github.com/ComparativeGenomicsToolkit/cactus/blob/master/doc/progressive.md)
it does mention running the alignment for like 2 months though so, intensive job... ("Both alignments were run on the AWS cloud over the course of 3 weeks for the avians and 2 months for the mammals")
they have also created an ecosystem for a data format called taffy that can store/index/decode multiple sequence alignments https://github.com/ComparativeGenomicsToolkit/taffy
See if these work:
https://sourceforge.net/projects/sam2fasta/files/
https://github.com/virus-evolution/gofasta (meant for COVID size genomes)