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
