Hello,

I have come across the same odd behaviour using mafsInRegion and from what I have noticed, the coordinates do not seem to matter, it's a problem with the number of lines in the bed file.

Running the program with a single line bed file returns the expected result. Running it with more than 5 lines will fetch the wrong region in some cases and will return an empty file for others.

Here's an example of a line that works when in a file by itself :

chr2    179456038       179462464       CATP00000883767.1       2       +

I get :

##maf version=1 scoring=autoMZ.v1
a score=0.000000
s hg19.chr2                        179456038 66 + 243199373 ATGGTAAGTACTGATGAGAAGTTGTCTGTTTCAACTGTGTAGTGCTCATCGGTTTTAATCTCACTC
s panTro2.chr2b                    183626670 66 + 248603653 ATGGTAAGTACTGATGAGAAGTTGTCTGTTTCAACTGTGTAGTGCTCATCGGTTTTAATCTCACTC
q panTro2.chr2b                                             999999999999999999999999999999999999999999999999999999999999999999
i panTro2.chr2b                    C 0 C 0
s gorGor1.Supercontig_0273903           2142 66 -      2212 
.....

which is the correct output.

when run in a file that contains 13 other lines and looks like this :

chr2    179406992       179407028       CATP00000004513.1       1       +
chr2    179436036       179440708       CATP00002698154.1       1       +
chr2    179442318       179442540       CATP00002345870.1       1       +
chr2    179464336       179464426       CATP00001176733.1       2       +
chr2    179500339       179500447       CATP00002315182.1       1       +
chr2    179474319       179474529       CATP00002157734.1       1       +
chr2    179629093       179638047       CATP00002900465.1       1       +
chr2    179479248       179479527       CATP00000210653.1       2       +
chr2    179436036       179462326       CATP00002016444.1       2       +
chr2    179439945       179440820       CATP00002345871.1       1       +
chr2    179404308       179404431       CATP00000332098.1       1       +
chr2    179482495       179482699       CATP00001007021.1       2       +
chr2    179412405       179412684       CATP00002348243.1       2       +
chr2    179456038       179462464       CATP00000883767.1       2       +

the coordinates get messed up and I get this :

 ##maf version=1 scoring=autoMZ.v1
a score=0.000000
s hg19.chr2                        179462327 96 + 243199373 AGACGCCTTGATGGCTCCTCTGGCAGTTCTTGATGACCATGG----ATGAGCTAATGGCTGTGGTCTCAATGGTGGCTTCTTGAGGTAAGGTTCTTTCAT
s panTro2.chr2b                    183632906 96 + 248603653 AGACGCCTTGATGGCTCCTCTGGCAGTTCTTGATGACCATGG----ATGAGCTAATGGCTGTGGTCTCAATGGTGGCTTCTTGAGGTAAGGTTCTTTCAT
q panTro2.chr2b                                             999999999999999999999999999999999999999999----999999999999999999999999999999999999999999999999999999
  ......

Notice how the start coordinates have changed in hg19.chr2 in the second case

When running the same multi-line file, the output file CATP00002345870.1.maf (corresponding to the bed file's first line) contains only :

##maf version=1 scoring=autoMZ.v1

However if I run the relevant line :

chr2    179406992       179407028       CATP00000004513.1       1       +

by itself, I get the correct output :

    ##maf version=1 scoring=autoMZ.v1
a score=0.000000
s hg19.chr2                        179406992 36 + 243199373 ATGCTGCTGCGACACTCTATGACCTCAGACTGCAAG
s anoCar1.scaffold_385                468376 36 +   1396552 ATTCTGCTTTTGCATTCTACGATTTCAGACTGTAGG
q anoCar1.scaffold_385                                      999999999999999999999999999999999999
i anoCar1.scaffold_385             I 8 C 0
s taeGut1.chr7                      17984185 36 +  39844632 ATTCTGCTCTTGCATTCTATAATTTCAGACTGCAGG
q taeGut1.chr7                                              999999999999999999999999999999999999
.....

I should specify that I am searching in the hg19/GRCh37 46-way MAF alignments downloaded from http://hgdownload.cse.ucsc.edu/goldenPath/hg19/multiz46way/

This feels like a bug to me, although I really hope it isn't. Any ideas?

Thanks

Hello, Thomas.

Thank you for questions about our utility "mafsInRegion".

Unfortunately, I was unable to replicate this issue. My input MAF file was hg38/multiz7way.maf and a BED file that contained the following two lines:

chr21   31659621        31668831        uc002ypa.4.1  
chr21   31659621        31668931        uc002ypa.4

As you can see both transcripts start at the same position, but have slightly different end positions. My output using this file was a file that contained MAF lines just as I would have expected:

##maf version=1 scoring=blastz
a score=0.000000
s hg38.chr21     31659621 547 +  46709983 GTTTGGGGCCAGAGTGGGCGAGGCGCGGAGGT-------CTGGCCTATAAAGTAGTCGCGGAGACGGGGTGCTG---GTTTGCGTCGTAGTCTCCTGCAGCG-------TCTGGGGTTTC------
...

Perhaps you can provide more information about what MAF files you are using with your example BED lines? I also attempted to replicate the issue using the BED lines that you provided and I didn't have any luck getting any meaningful output. I made a file with just a line for "transcript_1" and then I made another with both "transcript_1" and "transcript_2". In both cases, I used the hg38/multiz7way.maf MAF file I referenced before and I only got the output:

##maf version=1 scoring=blastz

Which would seem to indicate that there were no alignments in the MAF at this position? (I can also confirm this by looking at the Genome Browser for hg38 at this position with the 7-way multiz track active: https://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=mspeir&hgS_otherUserSessionName=hg38_biostars_226762.)

If you are still having issues, it would be helpful if you could post them to our Google Groups forum: https://groups.google.com/a/soe.ucsc.edu/forum/#!forum/genome-mirror, that way our whole team can see the question and help with an answer.

Thanks,

Matthew from the UCSC Genome Browser