Maf2SpanningSeq_PRANK.perl

[zuodabin]

Dear Professor:
The following error occurred when I was running Maf2SpanningSeq_PRANK.perl, I'm sure what went wrong, and how do I define elementID?

-- spanning coordinates after reading 8732282 blocks
Qspinosa is skipped because of inconsistencies (different chrom, strand, coordinate order)
Nanoranaparkeri is skipped because of inconsistencies (different chrom, strand, coordinate order)
wowa is skipped because of inconsistencies (different chrom, strand, coordinate order)
Ranatemporaria is skipped because of inconsistencies (different chrom, strand, coordinate order)

Using tempfile /dev/shm/prank.in.rb3scOlqi.fa ...
total number of bases to be aligned for HiC_scaffold_1: 0
only 0 species in the fasta file (/dev/shm/prank.in.rb3scOlqi.fa) --> nothing to align --> add string NOTHING_TO_ALIGN to /data/zuobin/comparative_genome/gemome_alignment/data/chr1BDB
*** ALL DONE ***

[MichaelHiller]

An alignment with Prank (or any other aligner) makes only sense if the sequence for this conserved element represents a single genomic locus in each species. The perl script extracts the coordinates of the cons element from a maf file. A maf file can have multiple maf blocks that make up the alignment of the conserved elements (with insertions in certain species in between maf blocks).
The perl script therefore checks whether the sequence coordinates in each maf block are co-linear in each species (meaning same strand, same order). If that is not the case (say 2 maf blocks and for wowa you have chr1 and chr8), then that species is removed, because the perl script does not know which sequence (chr1 or chr8??) to use and combining the two seqs makes no sense.

So in your case, you defined a super large region (entire scaffold instead of a conserved element??) and all species are removed because maf blocks are not co-linear. As a result, no species is left to align.

The likely solution is to run it over a set of conserved elements that typically represent a single locus.

[zuodabin]

Dear Professor,
I'm new to biological information, so I may not be able to understand what you mean, so I'll go withmafAddIRows->cat *.maf > >mafAddIRows.maf->mafIndex->mafExtract; After a process like mafExtract, which leads to an extract.out.maf file, how should I analyze the data next? I'm going to refer to an article published by your lab, **Convergent and lineage-specific genomic differences in limb regulatory elements in limbless reptile** Methods used in lineages. In addition, I have used Phastcons to get CNE.bed and exclude the code area, what should I do next?

Thank you professor

[MichaelHiller]

Pls run PhastCons or GERP on your genome-wide maf. That gives you conserved elements.
Our pipeline is meant to work with individual functional genomic regions, and conserved elements are a good proxy for that.

You probably fed the entire genome-wide maf into the perl script, which will never work.