Hi Michiel, One of our developers had this to say: "This may be an artifact of how we chop up the sequences for alignment on our cluster. rn4 chromosomes were chopped into 10,000,000-base chunks overlapping by 10,000 bases, and rheMac2 chroms were chopped into 50,000,000-base chunks with no overlap." Can you confirm that you also chopped up your sequences the same way?
It may also be helpful to you to display the differences between your file and the one in the browser. While we do not accept chain files as custom tracks, you can use the chainToPsl to make a PSL custom track for comparison. Please contact us again at [email protected] if you have any further questions. Best, Mary ------------------ Mary Goldman UCSC Bioinformatics Group On 7/4/11 10:20 PM, Michiel de Hoon wrote: > Many thanks for your help! With the extractRepeats script, I am getting > something closer to the alignment generated by UCSC, but there are still > discrepancies. > > These are the steps that I am taking to create the pairwise alignments, using > the alignment of rheMac2 chr7 to rn4 chr6 as an example: > > First I download chromFa.tar.gz and chromOut.tar.gz from > goldenPath/rheMac2/bigZips/ and goldenPath/rn4/bigZips/ > > For rheMac2, I run > > DateRepeats chr7.fa.out -query human -comp mouse -comp dog > > extractRepeats 1 chr7.fa.out_mus-musculus_canis-lupus-familiaris> > chr7.out.spec > > select_rpts chr7.out.spec 1 169801366 chr7.rpts > > strip_rpts chr7.fa chr7.rpts> chr7.strip.fa > > > For rn4, I run > > DateRepeats chr6.fa.out -query rat -comp human -comp mouse > > extractRepeats 1 chr6.fa.out_homo-sapiens_mus-musculus> chr6.out.spec > > select_rpts chr6.out.spec 1 147636619 chr6.rpts > > strip_rpts chr6.fa chr6.rpts> chr6.strip.fa > > > Then I run > > lastz chr7.strip.fa chr6.strip.fa --gap=400,30 --hspthresh=3000 > --gappedthresh=2200 --inner=2000 B=0> rheMac2.chr7.rn4.chr6.strip.plus.lav > > (using blastz instead of lastz gives the same result). > > Then I call restore_rpts, lavToPsl, and axtChain to generate > rheMac2.chr7.rn4.chr6.chain. > > This is the first chain in the rheMac2 chr7 to rn4 chr6 alignment that I > created: > > chain 526821446 chr7 169801366 + 87564296 168894851 chr6 147636619 + 64051113 > 138192842 1 > 17 2 0 > 57 0 1 > 19 0 1 > 7 7 0 > 139 6 0 > 15 1 0 > 12 7 0 > 4 3 0 > 11 0 1 > 57 1 0 > 52 0 1 > 23 6 0 > 61 0 1 > 51 2 0 > 71 1 0 > 14 1 0 > 15 4632 674 > 37 14 184<=== first discrepancy occurs here > ... > > and this is the corresponding chain generated by UCSC: > > chain 547645084 chr7 169801366 + 87564296 169142947 chr6 147636619 + 64051113 > 138454126 1 > 17 2 0 > 57 0 1 > 19 0 1 > 7 7 0 > 139 6 0 > 15 1 0 > 12 7 0 > 4 3 0 > 11 0 1 > 57 1 0 > 52 0 1 > 23 6 0 > 61 0 1 > 51 2 0 > 71 1 0 > 14 1 0 > 15 4632 674 > 39 0 154<=== first discrepancy occurs here > > This discrepancy can be identified already in > rheMac2.chr7.rn4.chr6.strip.plus.lav before calling restore_rpts, lavToPsl, > and axtChain, so it must be caused by some mistake in my call to lastz or > some earlier step. > > The only difference I can see with the UCSC pipeline is that I skip the call > to fasta-subseq, because I couldn't find this program and anyway I am running > lastz on the full chromosomes. > > Is there any other difference between the UCSC pipeline and what I am doing? > If not, would it be possible to make intermediate files of the UCSC pipeline > available for comparison? For example, it would be really informative if the > *.lav file with its header were available for one pairwise alignment for > comparison. > > Thanks again, and best wishes, > --Michiel de Hoon, RIKEN Omics Science Center. > > > --- On Fri, 6/24/11, Hiram Clawson<[email protected]> wrote: > >> From: Hiram Clawson<[email protected]> >> Subject: Re: [Genome] Repeat regions in whole-genome alignments >> To: "Michiel de Hoon"<[email protected]> >> Cc: [email protected] >> Date: Friday, June 24, 2011, 12:46 PM >> Good Morning Michiel: >> >> Please try this perl script: >> >> http://genomewiki.ucsc.edu/index.php/File:ExtractRepeats.txt >> >> --Hiram >> >> Michiel de Hoon wrote: >>> Dear Pauline, >>> >>> Many thanks for your reply. I think you are right and >> that I am inadvertently removing repeats that should not be >> removed. >>> To select the appropriate repeats from the DateRepeats >> output, UCSC uses the extractRepeats and extractLinSpecReps >> scripts. Can these scripts be made available somewhere? I >> couldn't find them in Jim Kent's software collection or >> elsewhere. >>> Thank you again again, >>> --Michiel > _______________________________________________ > Genome maillist - [email protected] > https://lists.soe.ucsc.edu/mailman/listinfo/genome _______________________________________________ Genome maillist - [email protected] https://lists.soe.ucsc.edu/mailman/listinfo/genome
