Maf file

[LemoAlex]

Dear Hiller Lab Team,

Thank you for sharing these amazing codes and resources !

I am trying to perform a Forward Genomics analysis on some fish dataset. I have a four species alignment maf file, which I obtained through progressive cactus (with a reference genome).

So, following the necessary steps, I first ran mafIndex:

mafIndex -chromSizes=~/reference.genome.chrom.sizes 4speciesAlignment.maf out.bb -keepBed=region.bdb

then mafExtract

mafExtract out.bb extract.out.maf -regionList=region.bdb -leaveEdgeMeta

These run without error and I obtain the demanded outputs.

However, when I then try :

perl Maf2SpanningSeq_PRANK.perl extract.out.maf Chromosome1* -runPrank -TwoBitSuffix ~/ForwardGenomics/genomes/ -treeFile tree.txt -BDBFile testBDB -noMissingSeq

• This is the first line of the reference.genome.chrom.size file

I get this error message:

there is no i-line that follows this s-line in the maf:
s Anc2.Anc2refChr14042 5866 234 - 7210 CACCAAAGACGGACAAGAAGAAACACAAAAAAGATCATCGACCAAACCAACTAACAAAAAGGGTATTAA [...]
s blenny.LR131933.1 22746305 234 - 30028995 CAGCCCTGATCGTGTACGCCCACCAGAGTCCGGGTTCGTTCAATGCAGCGGTGCGTGACGTAGCGATGCAG [...]

Do you have any idea of why this is not working ? Could it be due to the fact that I am using a maf file which was not obtained through Lastz? I compared it and the structure look the same so I am a bit confused about this.

Thank you for your time and help,

Best Regards,

Alexandre

[MichaelHiller]

Hi Alexandre,

pls run mafAddIRows on the maf to get the i and e lines. Then mafIndex etc. Then it will work.

Michael

[LemoAlex]

Dear Prof Hiller,

Thank you for your answer.

I was able to run mafAddIRows, but it required me to run mafSplit before. After that, I obtained one maf file per region (so thousands of Mafs), for which I have to run mafAddIRows individually. This seems odd, and I am not sure if it's the right way to do it, since normally I should run mafIndex & mafExtract on the big alignment file ?

Am I mistaken or is something that I am not getting straight?

Thank you for your help,

Alexandre

[MichaelHiller]

Good point. We are running mafAddIRows already on the split maf files
You can use mafSplit -byTarget (maybe with -useSequenceName, I don't remember).
Then run mafAddIRows. Then combine again (strip the header to have this only once) and then mafIndex.
You can also index the per-chrom mafs of course.

Michael

[LemoAlex]

Thank you for your answer.

So I ran mafSplit and then maffAddI rows, for each region (scaffold or chromosome) of the reference genome, and afterwards simply cat all the files. Now I should be able to run Maf2SpanningSeq_PRANK.perl.

However, when running maffAddIrows, I got an error message for the majority regions (i.e. ~80% of the maf files) :

zero left status species.region:0

and obtain empty files.

I thought it was possible that some regions are not indexable (?), but even some big regionfor which the initial .maf file is rather big get me this message. I am confused of what is happening here

Thanks again,
Best Regards,
Alexandre

[MichaelHiller]

I haven't seen this error. Does the maf file look OK? Why do some maf files work and others don't?

[LemoAlex]

The mafs look ok to me . I don't know why some work and others don't.

Here is the beggining of a working file :

##maf version=1 scoring=blastz
a score=0.000000
s blenny.CAAAFJ010000002.1 0 16 + 13669 ANNNNNNNNNNNNNNN
s blenny.CAAAFJ010000002.1 320 10 + 13669 TGTTTGGTGT
s Anc2.Anc2refChr3323465 1 2 - 3 TC--------
s Anc2.Anc2refChr21310 642 8 + 2563 --TTTTGGTT
s stickleback.AANH01000931.1 47440 10 + 93272 TGTTTGGTGT

a score=0.000000
s blenny.CAAAFJ010000002.1 330 20 + 13669 TATTTATACACACAAATTAA
s Anc2.Anc2refChr21310 650 20 + 2563 GGTTGTTCTTGGGCTCTCGT
s stickleback.AANH01000931.1 47451 20 + 93272 TATTCAAACACACAACTGAA

and the beggining of a non-working file:

##maf version=1 scoring=blastz
a score=0.000000
s blenny.CAAAFJ010000001.1 0 9 + 37463 CATAAACAA
s Anc2.Anc2refChr12290 173 9 - 312 ACACATGTA
s blenny.CAAAFJ010000237.1 102605 9 - 111081 CATAAACAA
s blenny.CAAAFJ010000237.1 32179 9 - 111081 CATAAACAA
s blenny.CAAAFJ010000001.1 1330 9 - 37463 CATAAACAA
s blenny.LR131939.1 8120627 9 - 30072547 CATAAACAA

a score=0.000000
s blenny.CAAAFJ010000001.1 1626 29 + 37463 TTTTTGTGACACCAACCCTCTTGTGCTCC
s Anc2.Anc2refChr8970 331 29 - 1392 GACGAGCGATACAAAGCAATCCAAAAACA
s blenny.CAAAFJ010000237.1 104236 29 - 111081 TTTTTGTGACACCAACCCTCTTGTGCTCC
s blenny.CAAAFJ010000237.1 33801 29 - 111081 TTTTTGTGACACCAACCCTCTTGTGCTCC
s blenny.LR131939.1 8122236 29 - 30072547 TTTTTGTGACACCAACCCTCTTGTGCTCC
s stickleback.AANH01014576.1 2703 29 - 4163 GTCATGTGACACCGACCCTCCCGTCCCCC

I can't see any real difference there, but maybe I am missing something. The reference genome has 322 different regions ( so 322 files obtained from mafSplit), and there does not seem to be a pattern for the failing mafs (as there are both big and small files/regions which are failing)

Thank you for your time,
Alexandre

[MichaelHiller]

Do you have a 2bit file for each of the species? I think this is what mafAddIRows needs.
Also, is this a proper multiz maf file? Anc sounds like reconstructed ancestors, which multiz would not do.

If you do have 2bit files and it keeps crashing, pls send the error message. But not sure I can debug this, as mafAddIRows is a UCSC tool (they have a forum as well)

[LemoAlex]

Hello again,

I finally found out where the issue came from, or the first part at least. Maybe this can be useful to someone in the future. Basically, my alignment was obtained through cactus and not through lastz/multiz pipeline. Therefore, in the hal2maf conversion, I had to add the -noAncestor and -noDupes (remove paralog edges) option for the maf to be in a suitable format, i.e.:

##maf version=1 scoring=N/A
a
s spec.LR131921.1 1000 1000 + 25704503
tatgcatttatgtaaatgtat[...]
a
s spec.LR131921.1 2000 1000 + 25704503
gtgtaatacacgtgcgctt[...]

and so on...

So I am able to run mafSplit -> mafAddIrows --> cat all the .maf files --> mafIndex -->mafExtract, using the regionList obtained in mafIndex (keepBed).

Then, when I am trying to run Maf2SpanningSeq_PRANK.perl, I get a new error (I added the -verbose, to get more info).

##maf version=1 scoring=mafExtract
a score=0.000000
s ref.species.LR131921.1 9726354 558 + 25704503 TGAAATGTAAATTGTTAACAATATAAGCTACACACTCACTTCTCCAGCTAAAGTGAATTTGAACTTTGTAAAAGGTTAGGAAACAGTTCAACACGTCATAAATCAGTGGAAGTAGAAATAGCTGATGAAACGTATGTCACAATGTGTCTCACAGAAACATGATCATAATTAATAACACGTTAAATTAAAAAGGTTTATGAACTCATTTATAAAAGAAAAAGCTTAATATGCAGCATGGTGCTGCATCATTAAACCATTGCACAACACCTTCAGTATAGTTCCTATTTAGATAGGATATAACACACGCTGTAAAACATTATGAAGCCATCAGAACAATAGTGCATCCGTTTCCTGACCTAAAAATAATATTTAAAAGTTCAAGCTCATATGCTCTCAATGCTACCGATATGAATCCTGTTTCATGTTTGAGTTTTTACTTGTTAACACATCTGATGTCACTTGTGGGTGCAGTATTCTATAACAATGGCAGTACAGTAAATAAAGTGGGAACCTGCTTTAGAAATAGAAAATGTCCCCACACTGACCCCGAACAGGCTG
e species2.AANH01000521.1 130795 1 + 207688 n
--> reference species #mafblocks: 1
-- spanning coordinates after reading 1 blocks
ref.species LR131921.1:9726354-9726912 + srcSize 25704503 span 558
species2 AANH01000521.1:130795-130796 + srcSize 207688 span 1 !!!! ONLY E-line blocks given --> seq will not be in the output !!!
sh: 1: set: Illegal option -o pipefail
ERROR in getSeqTwoBitFile: set -o pipefail; echo -e "LR131921.1\t9726354\t9726912\tname\t1000\t+" | twoBitToFa /home/alexandre/ForwardGenomics/genomes/refspec/refspec.2bit stdout -bed=stdin -noMask failed

Any idea of what is going wrong? If I simply run

echo -e "LR131921.1\t9726354\t9726912\tname\t1000\t+" | twoBitToFa /home/alexandre/ForwardGenomics/genomes/refspec/refspec.2bit stdout -bed=stdin -noMask

Then I get this output:

name
TGAAATGTAAATTGTTAACAATATAAGCTACACACTCACTTCTCCAGCTA
AAGTGAATTTGAACTTTGTAAAAGGTTAGGAAACAGTTCAACACGTCATA
AATCAGTGGAAGTAGAAATAGCTGATGAAACGTATGTCACAATGTGTCTC
ACAGAAACATGATCATAATTAATAACACGTTAAATTAAAAAGGTTTATGA
ACTCATTTATAAAAGAAAAAGCTTAATATGCAGCATGGTGCTGCATCATT
AAACCATTGCACAACACCTTCAGTATAGTTCCTATTTAGATAGGATATAA
CACACGC

So I don't know what is happening since this command seem to be working.

I am really sorry to bother you again with all these issues!

Thank you for your help,

Alexandre

[MichaelHiller]

Hi Alex,
Yes, I was suspecting that this is a cactus alignment.
If I understand correctly, the issue is with pipefail.
This was simply a safety measure for bash to assure a proper error exit if any command in the pipe (e.g. twoBitToFa) fails.

If that doesn't work for you (using csh??), then pls remove the "set -o pipefail" parts from the perl script in all commands.
Hope it works then.

Best
Michael

[LemoAlex]

Dear Michael,

Thank you for your input once again. I removed the set-o pipefail from the script. However, i don't know what is happening, but Maf2SpanningSeq is still failing:

perl Maf2SpanningSeq_PRANK.perl 2182.maf 2182 -treeFile tree.nh -BDBFile chr1BDB -noMissingSeq -verbose

##maf version=1 scoring=mafExtract
a score=0.000000
s refSpec.LR131916          1000 684 + 27055867 
TAACAGAGGAATGCCAAACAGGCG
s stickleback.AANH01012136  868 112 -    64769 
GTAACACTGACCCTCTGAACTCTGAC
i stickleback.AANH01012136 I 125 I 3
--> reference species blenny   #mafblocks: 1
-- spanning coordinates after reading 1 blocks
              blenny    LR131916:1000-1684 +            srcSize 27055867                span 684
         stickleback    AANH01012136:868-980 -          srcSize 64769           span 112                +coords: AANH01012136:63789-63901
>
invalid unsigned integer: "LR131916"
ERROR in getSeqTwoBitFile: echo -e "LR131916\t1000\t1684\tname\t1000\t+" | twoBitToFa /home/alexandre/ForwardGenomics/genomes/blenny/blenny.2bit stdout -bed=stdin -noMask failed

What is this "unsigned integer" error exactly? I checked and the chromosome names are totally matching with the files obtained from twobitinfo (i.e. corresponding to the chrom.size file).

I also tried to compare all the bed /maf /tree files with the ones from your example files, and the structure is matching, so the input files should be just fine.

Sorry for the bother and thank you for your help,

Best,
Alexandre

[MichaelHiller]

Pls try this
echo -e "LR131916\t1000\t1684\tname\t1000\t+" > test.bed
twoBitToFa /home/alexandre/ForwardGenomics/genomes/blenny/blenny.2bit stdout -bed=test.bed

If that still fails, pls do
twoBitInfo /home/alexandre/ForwardGenomics/genomes/blenny/blenny.2bit chrom.sizes
and check whether LR131916 occurs in the genome (chrom.sizes list).

[LemoAlex]

On their own, these work :

echo -e "LR131916\t1000\t1684\tname\t1000\t+" > test.bed
twoBitToFa /home/alexandre/ForwardGenomics/genomes/blenny/blenny.2bit stdout -bed=test.bed

name
[...]ACACAGACGCTCTGCATGGATGTGATCTGACGGAAGACTTTATCAGCAGCA

I tried replacing the original code in the MyKentFunctions.pm file:

sub getSeqTwoBitFile {
	my ($chr, $start, $end, $strand, $TwoBitFile, $unmask) = @_;

	# run twoBitToFa with -noMask if unmask set to 1
	my $unmaskParameter = "";
	$unmaskParameter = "-noMask" if ($unmask == 1);

	if (! -e $TwoBitFile) {
		die "ERROR in getSeqTwoBitFile: $TwoBitFile does not exist\n";
	}
	
	my $call = "echo -e \"$chr\\t$start\\t$end\\tname\\t1000\\t$strand\" > test.bed
		 twoBitToFa $TwoBitFile stdout -bed=test.bed $unmaskParameter";
	print "$call\n" if ($verbose);
	my @res = `$call`;
	die "ERROR in getSeqTwoBitFile: $call failed\n" if ($? != 0 || ${^CHILD_ERROR_NATIVE} != 0);
	print STDERR "ERROR: cannot get seq by calling: $call\n" if ($#res < 1);
	
	# skip first fasta header and remove the \n between the lines
	shift @res; $" = ""; chomp @res; my $seq = "@res";

	# if unmask set to 2, replace lower case letters by N
	if ($unmask == 2) {
		$seq =~ s/[acgtn]/N/g;
	}
	return $seq;
}

But I still get the same error:

invalid unsigned integer: "LR131916"
ERROR in getSeqTwoBitFile: echo -e "LR131916\t1000\t1684\tname\t1000\t+" > test.bed
                 twoBitToFa /home/alexandre/ForwardGenomics/genomes/blenny/blenny.2bit stdout -bed=test.bed -noMask failed

I also checked and LR131916 occurs in the chrom sizes list.

I am sorry to bother you with this problem, but I am really running out of ideas on what is happening here.

Best,
Alexandre

[MichaelHiller]

Hmm, I am also running out of ideas. The code was worked millions of times and if you can manually run the command that fails on the shell, then I don't know why $call fails.
Is this the only genomic region that crashes? What if you run it on a region for another reference chrom (not LR131916)?

Can you try to outcomment the "die "ERROR" command which causes the exit and add a
print "==> RESULT: @res\n" to see whether the command actually fetches the right sequence from the 2bit file?

[LemoAlex]

Hi,

Good news, I managed to resolve the issue. I don't exactly know why it was failing, but I thought the problem could come from the echo -e function. I replaced it with "printf" as in :

`my $call="printf \"$chr\\t$start\\t$end\\tname\\t1000\\t$strand\"|twoBitToFa $TwoBitFile stdout -bed=stdin $unmaskParameter";	print "$call\n" if ($verbose);`

And now I am able to run Maf2SpanningSeq until the end, where I get the message

> 
Write result to testBDB with the key: 2182 ... DONE
*** ALL DONE ***

I will now try the rest of the pipeline :)

Thank you for your help!

-Alexandre

[MichaelHiller]

Wow, I have no idea why printf but not echo works.
Glad this is resolved