# for emacs: -*- mode: sh; -*- # $Id: melGal1.txt,v 1.6 2010/05/06 16:27:44 chinhli Exp $ # Meleagris gallopavo (turkey) -- Turkey Genome Consortium Turkey_2.01 (2009-12-10) ATE: 10-DEC-2009 ORGANISM: Meleagris gallopavo TAXID: 9103 ASSEMBLY LONG NAME: Turkey_2.01 ASSEMBLY SHORT NAME: Turkey_2.01 ASSEMBLY SUBMITTER: Turkey Genome Consortium ASSEMBLY TYPE: Haploid NUMBER OF ASSEMBLY-UNITS: 1 Assembly Accession: GCA_000146605.1 # file template copied from susScr2.txt # Meleagris gallopavo (NCBI Project ID: 10805, Accession: GCA_000146605.1) # by Turkey Genome Consortium # sequence: # ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/Meleagris_gallopavo/Turkey_2.01/ # http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=ADDD00 ########################################################################## # Download sequence (DONE - 2010-10-22 - Chin) mkdir /hive/data/genomes/melGal1 cd /hive/data/genomes/melGal1 mkdir genbank cd genbank wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \ --no-remove-listing -np \ "ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/Meleagris_gallopavo/Turkey_2.01/*" # FINISHED --2010-10-22 08:28:53-- # Downloaded: 243 files, 611M in 9m 19s (1.09 MB/s) # Read ASSEMBLY_INFO mkdir ucscChr # stay at genbank directory # fixup the accession names to become UCSC chrom names export S=Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` echo "${ACC} -> chr${C}" zcat ${S}/AGP/chr${C}.comp.agp.gz \ | sed -e "s/^${ACC}/chr${C}/" | gzip > ucscChr/chr${C}.agp.gz done export S=Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` echo "${ACC} -> chr${C}" echo ">chr${C}" > ucscChr/chr${C}.fa zcat ${S}/FASTA/chr${C}.fa.gz | grep -v "^>" >> ucscChr/chr${C}.fa gzip ucscChr/chr${C}.fa & done # Check them with faSize faSize Primary_Assembly/assembled_chromosomes/FASTA/chr*.fa.gz # 2784748484 bases (1600136831 N's 1184611653 real 1184611653 upper # 0 lower) in 27 sequences in 27 files faSize ucscChr/chr*.fa.gz # 1061800384 bases (125894715 N's 935905669 real 935905669 upper # 0 lower) in 5890 sequences in 33 files # For unplaced scalfolds, named them as chrUn_xxxxxxxx # where xxxxxx is the original access id as: chrUn_GL340781.1 # The ".1" at the end need to be filter out since # MySQL does not allow "." as part of the table name and # will casue problems in genbank task step later export S=Primary_Assembly/unplaced_scaffolds zcat ${S}/AGP/unplaced.scaf.agp.gz | grep "^#" > ucscChr/chrUn.agp # append the gap records zcat ${S}/AGP/unplaced.scaf.agp.gz | grep -v "^#" \ | sed -e "s/^/chrUn_/" -e "s/\.1//" >> ucscChr/chrUn.agp gzip ucscChr/chrUn.agp & zcat ${S}/FASTA/unplaced.scaf.fa.gz \ | sed -e "s#^>.*|gb|#>chrUn_#; s#|.*##" -e "s/\.1//" \ | gzip > ucscChr/chrUn.fa.gz # about 5858 sequences in the unplaced zcat ucscChr/chrUn.fa.gz | grep "^>" | wc # 5858 5858 93728 # Check them with faSize faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz # 21513314 bases (2885592 N's 18627722 real 18627722 upper # 0 lower) in 5858 sequences in 1 files # 75747883 bases (59104875 N's 16643008 real 16643008 upper 0 lower) # in 1190 sequences in 1 files faSize ucscChr/chrUn.fa.gz # 21513314 bases (2885592 N's 18627722 real 18627722 upper # 0 lower) in 5858 sequences in 1 files # N50 mkdir N50 export S=Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` echo ">chr${C}" > ucscChr/chr${C}.fa zcat ucscChr/chr${C}.fa.gz >> N50/melGal1.fa done zcat ucscChr/chrUn.fa.gz >> N50/melGal1.fa faCount N50/melGal1.fa | awk ' /^chr/ {print $1, $2}' > N50/chrom.sizes n50.pl N50/chrom.sizes # reading: N50/chrom.sizes # contig count: 5890, total size: 1061800384, one half size: 530900192 # cumulative N50 count contig contig size 502447651 4 chrZ 81011772 530900192 one half size 577312103 5 chr4 74864452 rm N50/melGal1.fa ######################################################################### # Initial makeGenomeDb.pl (DONE - 2010-11-03 - Chin) cd /hive/data/genomes/melGal1 cat << '_EOF_' > melGal1.config.ra # Config parameters for makeGenomeDb.pl: db melGal1 clade vertebrate genomeCladePriority 55 scientificName Meleagris gallopavo commonName Turkey assemblyDate Dec. 2009 assemblyLabel TGC (NCBI Project ID: 10805, Accession: GCA_000146605.1) assemblyShortLabel TGC Turkey_2.01 orderKey 436 mitoAcc NC_010195 fastaFiles /hive/data/genomes/melGal1/genbank/ucscChr/chr*.fa.gz agpFiles /hive/data/genomes/melGal1/genbank/ucscChr/chr*.agp.gz # qualFiles none dbDbSpeciesDir turkey taxId 9103 '_EOF_' # << happy emacs time makeGenomeDb.pl -noGoldGapSplit -workhorse=hgwdev melGal1.config.ra \ > makeGenomeDb.log 2>&1 & # real 7m58.244s # add the trackDb entries to the source tree, and the 2bit link: ln -s `pwd`/melGal1.unmasked.2bit /gbdb/melGal1/melGal1.2bit # Per instructions in makeGenomeDb.log: # cd ~/kent/src/hg/makeDb/trackDb # edit makefile to add melGal1 to DBS. # git add turkey/melGal1/*.{ra,html} # git commit -m "Added melGal1 to DBS." makefile # git commit -m "Initial descriptions for melGal1." # turkey/melGal1/*.{ra,html} # git pull; git push # Run make update DBS=melGal1 and make alpha when done. # (optional) Clean up /cluster/data/melGal1/TemporaryTrackDbCheckout ######################################################################### # RepeatMasker (DONE - 2010-11-04 - Chin) mkdir /hive/data/genomes/melGal1/bed/repeatMasker cd /hive/data/genomes/melGal1/bed/repeatMasker time nice -n +19 doRepeatMasker.pl -buildDir=`pwd` \ -workhorse=hgwdev -bigClusterHub=swarm -noSplit melGal1 > do.log 2>&1 & # real 36m27.415s cat faSize.rmsk.txt # 1061817101 bases (125894715 N's 935922386 real 874501528 upper # 61420858 lower) in 5891 sequences in 1 files # %5.78 masked total, %6.56 masked real ######################################################################### # simpleRepeats ( DONE - 2010-11-05 - Chin) mkdir /hive/data/genomes/melGal1/bed/simpleRepeat cd /hive/data/genomes/melGal1/bed/simpleRepeat time nice -n +19 doSimpleRepeat.pl -buildDir=`pwd` -workhorse=hgwdev \ -bigClusterHub=swarm -smallClusterHub=swarm melGal1 > do.log 2>&1 & # real 7m12.044s cat fb.simpleRepeat # 4997546 bases of 935922386 (0.534%) in intersection # add to the repeatMasker cd /hive/data/genomes/melGal1 twoBitMask melGal1.rmsk.2bit -add bed/simpleRepeat/trfMask.bed melGal1.2bit # safe to ignore warnings about >=13 fields twoBitToFa melGal1.2bit stdout | faSize stdin > melGal1.2bit.faSize.txt cat melGal1.2bit.faSize.txt # 1061817101 bases (125894715 N's 935922386 real 874132579 upper # 61789807 lower) in 5891 sequences in 1 files # %5.82 masked total, %6.60 masked real ######################################################################### # Marking *all* gaps - they are all in the AGP file # (DONE 2010-11-05 - Chin) mkdir /hive/data/genomes/melGal1/bed/allGaps cd /hive/data/genomes/melGal1/bed/allGaps time nice -n +19 findMotif -motif=gattaca -verbose=4 \ -strand=+ ../../melGal1.unmasked.2bit > findMotif.txt 2>&1 # real 0m22.247s grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed featureBits melGal1 -not gap -bed=notGap.bed # 935922386 bases of 935922386 (100.000%) in intersection featureBits melGal1 allGaps.bed notGap.bed -bed=new.gaps.bed # 0 bases of 935922386 (0.000%) in intersection # what is the highest index in the existing gap table: hgsql -N -e "select ix from gap;" melGal1 | sort -n | tail -1 # 53110 hgsql -e "select count(*) from gap;" melGal1 # 146748 ######################################################################## # Create kluster run files (DONE 2010-11-05 - Chin) # numerator is melGal1 gapless bases "real" as reported by: featureBits -noRandom -noHap melGal1 gap # 123009123 bases of 917294664 (13.410%) in intersection # 1600136831 bases of 1184628269 (135.075%) in interddsection # denominator is hg19 gapless bases as reported by: # featureBits -noRandom -noHap hg19 gap # 234344806 bases of 2861349177 (8.190%) in intersection # 1024 is threshold used for human -repMatch: calc \( 917294664 / 2861349177 \) \* 1024 # ( 917294664 / 2861349177 ) * 1024 = 328.275117 # ==> use -repMatch=400 according to size scaled down from 1024 for human. # and rounded down to nearest 50 cd /hive/data/genomes/melGal1 blat melGal1.2bit \ /dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/melGal1.11.ooc \ -repMatch=300 & # Wrote 22129 overused 11-mers to jkStuff/melGal1.11.ooc mkdir /hive/data/staging/data/melGal1 cp -p melGal1.2bit jkStuff/melGal1.11.ooc /hive/data/staging/data/melGal1 cp -p chrom.sizes /hive/data/staging/data/melGal1 # check non-bridged gaps to see what the typical size is: hgsql -N \ -e 'select * from gap where bridge="no" order by size;' melGal1 \ | sort -k7,7nr # most gaps have size > 100,000 # decide on a minimum gap for this break gapToLift -verbose=2 -minGap=20000 melGal1 jkStuff/nonBridged.lft \ -bedFile=jkStuff/nonBridged.bed cp -p jkStuff/nonBridged.lft \ /hive/data/staging/data/melGal1/melGal1.nonBridged.lft # ask cluster-admin to copy (evry time if any file chsnged) # /hive/data/staging/data/melGal1 directory to cluster nodes # /scratch/data/melGal1 ######################################################################## # GENBANK AUTO UPDATE (DONE - 201w1-02-16 - Chin) # 2011-10-07 Note: Re-configure per Mark's comments: # genbank should not have native refseq enabled for this genome, as there are # no native refseqs. Xeno refseqs should be enabled. Please fix configuration # file and do a drop and reload. ssh hgwdev cd $HOME/kent/src/hg/makeDb/genbank git pull # edit etc/genbank.conf to add melGal1 just after chicken galGal2 # melGal1 (turkey) melGal1.serverGenome = /hive/data/genomes/melGal1/melGal1.2bit melGal1.clusterGenome = /scratch/data/melGal1/melGal1.2bit melGal1.ooc = /scratch/data/melGal1/melGal1.11.ooc melGal1.lift = no melGal1.perChromTables = no melGal1.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter} melGal1.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter} melGal1.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter} melGal1.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter} melGal1.genbank.est.xeno.pslCDnaFilter = ${ordered.genbank.est.xeno.pslCDnaFilter} melGal1.downloadDir = melGal1 melGal1.refseq.mrna.native.load = no melGal1.refseq.mrna.xeno.load = yes melGal1.refseq.mrna.xeno.loadDesc = yes melGal1.genbank.mrna.native.load = yes melGal1.genbank.mrna.native.loadDesc = yes melGal1.genbank.mrna.xeno.load = yes melGal1.genbank.mrna.xeno.loadDesc = yes melGal1.genbank.est.native.load = yes melGal1.genbank.est.native.loadDesc = yes git add etc/genbank.conf git commit -m "Added melGal1" etc/genbank.conf git pull git push # update /cluster/data/genbank/: make etc-update # Edit src/lib/gbGenome.c to add new species. With these two lines: # static char *melGalNames[] = {"Meleagris gallopavo", NULL}; # ... later ... # {"melGal", melGal1Names}, # gbGenome.c is in # /cluster/home/chinhli/kent/src/hg/makeDb/genbank/src/lib # make and checkin make install-server git add src/lib/gbGenome.c git commit -m "adding melGal1 Turkey" src/lib/gbGenome.c git pull git push ssh genbank screen # control this business with a screen since it takes a while cd /cluster/data/genbank time nice -n +19 ./bin/gbAlignStep -initial melGal1 & # logFile: var/build/logs/2011.02.16-09:21:59.melGal1.initalign.log # real 109m5.431s # load database when finished ssh hgwdev cd /cluster/data/genbank time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad melGal1 & # logFile: var/dbload/hgwdev/logs/2011.02.16-12:45:15.dbload.log # real 23m10.935s # Re-run load step 2011-10-20: # logFile: var/dbload/hgwdev/logs/2011.10.20-10:10:19.dbload.log # real 53m49.262s # enable daily alignment and update of hgwdev cd ~/kent/src/hg/makeDb/genbank git pull # add melGal1 to: etc/align.dbs etc/hgwdev.dbs git add etc/align.dbs git add etc/hgwdev.dbs git commit -m "Added melGal1 - Turkey" etc/align.dbs etc/hgwdev.dbs git pull git push make etc-update ######################################################################### # reset default position as galGal3 Chicken (DONE 2011-03-08 - Chin) # hgsql -e \ 'update dbDb set defaultPos="chr5:57710001-57780000" where name="melGal1";' \ hgcentraltest ############################################################################ # ctgPos2 track - showing clone sequence locations on chromosomes # (DONE 2011-03-08 - Chin) # NOTE - create melGal1 entry in all.joiner since this is a new species mkdir /hive/data/genomes/melGal1/bed/ctgPos2 cd /hive/data/genomes/melGal1/bed/ctgPos2 cat << '_EOF_' > agpToCtgPos2.pl #!/usr/bin/env perl use warnings; use strict; my $argc = scalar(@ARGV); if ($argc != 1) { printf STDERR "usage: zcat your.files.agp.gz | agpToCtgPos2.pl /dev/stdin > ctgPos2.tab\n"; exit 255; } my $agpFile = shift; open (FH, "<$agpFile") or die "can not read $agpFile"; while (my $line = ) { next if ($line =~ m/^#/); chomp $line; my @a = split('\s+', $line); next if ($a[4] =~ m/^N$/); next if ($a[4] =~ m/^U$/); my $chrSize = $a[2]-$a[1]+1; my $ctgSize = $a[7]-$a[6]+1; die "sizes differ $chrSize != $ctgSize\n$line\n" if ($chrSize != $ctgSize); printf "%s\t%d\t%s\t%d\t%d\t%s\n", $a[5], $chrSize, $a[0], $a[1]-1, $a[2], $a[4]; } close (FH); '_EOF_' # << happy emacs chmod +x agpToCtgPos2.pl export S=../../genbank/Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` zcat ${S}/AGP/chr${C}.agp.gz \ | sed -e "s/^${ACC}/chr${C}/" done | ./agpToCtgPos2.pl /dev/stdin > ctgPos2.tab hgLoadSqlTab melGal1 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab ############################################################################ # running cpgIsland business (DONE -2011-03-08 - Chin) mkdir /hive/data/genomes/melGal1/bed/cpgIsland cd /hive/data/genomes/melGal1/bed/cpgIsland cvs -d /projects/compbio/cvsroot checkout -P hg3rdParty/cpgIslands cd hg3rdParty/cpgIslands # needed to fixup this source, adding include to readseq.c: #include "string.h" # and to cpg_lh.c: #include "unistd.h" #include "stdlib.h" # and fixing a declaration in cpg_lh.c sed -e "s#\(extern char\* malloc\)#// \1#" cpg_lh.c > tmp.c mv tmp.c cpg_lh.c make cd ../../ ln -s hg3rdParty/cpgIslands/cpglh.exe mkdir -p hardMaskedFa cut -f1 ../../chrom.sizes | while read C do echo ${C} twoBitToFa ../../melGal1.2bit:$C stdout \ | maskOutFa stdin hard hardMaskedFa/${C}.fa done ssh swarm cd /hive/data/genomes/melGal1/bed/cpgIsland mkdir results cut -f1 ../../chrom.sizes > chr.list cat << '_EOF_' > template #LOOP ./runOne $(root1) {check out exists results/$(root1).cpg} #ENDLOOP '_EOF_' # << happy emacs # the faCount business is to make sure there is enough sequence to # work with in the fasta. cpglh.exe does not like files with too many # N's - it gets stuck cat << '_EOF_' > runOne #!/bin/csh -fe set C = `faCount hardMaskedFa/$1.fa | grep ^chr | awk '{print $2 - $7 }'` if ( $C > 200 ) then ./cpglh.exe hardMaskedFa/$1.fa > /scratch/tmp/$1.$$ mv /scratch/tmp/$1.$$ $2 else touch $2 endif '_EOF_' # << happy emacs chmod 775 runOne gensub2 chr.list single template jobList para create jobList para try para check ... etc para time para problems para status # then, kick it with para push # check it with plb # when all are done, para time shows: # para time # 5891 jobs in batch # 130198 jobs (including everybody's) in Parasol queue or running. # Checking finished jobs # Completed: 5891 of 5891 jobs # CPU time in finished jobs: 86s 1.44m 0.02h 0.00d 0.000 y # IO & Wait Time: 19261s 321.01m 5.35h 0.22d 0.001 y # Average job time: 3s 0.05m 0.00h 0.00d # Longest finished job: 26s 0.43m 0.01h 0.00d # Submission to last job: 369s 6.15m 0.10h 0.00d # Transform cpglh output to bed + catDir results | awk '{ $2 = $2 - 1; width = $3 - $2; printf("%s\t%d\t%s\t%s %s\t%s\t%s\t%0.0f\t%0.1f\t%s\t%s\n", $1, $2, $3, $5,$6, width, $6, width*$7*0.01, 100.0*2*$6/width, $7, $9); }' > cpgIsland.bed ssh hgwdev cd /hive/data/genomes/melGal1/bed/cpgIsland hgLoadBed melGal1 cpgIslandExt -tab \ -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed # Loaded 13372 elements of size 10 # Sorted # Creating table definition for cpgIslandExt # Saving bed.tab # Loading melGal1 # cleanup rm -fr hardMaskedFa ######################################################################### # all.joiner update, downloads and in pushQ - ( DONE 2011-03-08 - Chin) cd $HOME/kent/src/hg/makeDb/schema # fixup all.joiner until this is a clean output # add melGal after galGal joinerCheck -database=melGal1 -all all.joiner mkdir /hive/data/genomes/melGal1/goldenPath cd /hive/data/genomes/melGal1/goldenPath makeDownloads.pl melGal1 > do.log 2>&1 # now ready for pushQ entry mkdir /hive/data/genomes/melGal1/pushQ cd /hive/data/genomes/melGal1/pushQ makePushQSql.pl melGal1 > melGal1.pushQ.sql 2> stderr.out # check for errors in stderr.out, some are OK, e.g.: # WARNING: hgwdev does not have /gbdb/melGal1/wib/gc5Base.wib # WARNING: hgwdev does not have /gbdb/melGal1/wib/quality.wib # WARNING: hgwdev does not have /gbdb/melGal1/bbi/quality.bw # WARNING: melGal1 does not have seq # WARNING: melGal1 does not have extFile # # WARNING: Could not tell (from trackDb, all.joiner and hardcoded lists of # supporting and genbank tables) which tracks to assign these tables to: # ctgPos2 # copy it to hgwbeta scp -p melGal1.pushQ.sql hgwbeta:/tmp ssh hgwbeta cd /tmp hgsql qapushq < melGal1.pushQ.sql # in that pushQ entry walk through each entry and see if the # sizes will set properly ########################################################################## # BLATSERVERS ENTRY (DONE 2011-03-23 - Chin) # After getting a blat server assigned by the Blat Server Gods, # # blat13 # trans port 17828 # untrans port 17829 hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("melGal1", "blat13", "17828", "1", "0"); \ INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("melGal1", "blat13", "17829", "0", "1");' \ hgcentraltest # test it with some sequence ######################################################################### # SWAP melGal1 lastz (DONE - 2011-03-24 - Chin) # original alignment cat fb.galGal3.chainMelGal1Link.txt # 924755061 bases of 1042591351 (88.698%) in intersection # and the swap mkdir /hive/data/genomes/melGal1/bed/blastz.galGal3.swap cd /hive/data/genomes/melGal1/bed/blastz.galGal3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/galGal3/bed/lastzMelGal1.2011-03-24/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 38m4.406s cat fb.melGal1.chainGalGal3Link.txt # 874579772 bases of 935922386 (93.446%) in intersection ######################################################################### # BLASTZ/CHAIN/NET Turkey vs Zebra Fish taeGut1 (DONE - 2011-03-28 - Chin) mkdir /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28 cd /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28 cat << '_EOF_' > DEF # Zebra Finch vs Turkey # TARGET: Turkey melGal1 SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit SEQ1_LEN=/scratch/data/melGal1/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 # QUERY: Zebra Finch taeGut1 - single chunk big enough to run entire chrom SEQ2_DIR=/scratch/data/taeGut1/taeGut1.2bit SEQ2_LEN=/scratch/data/taeGut1/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 BASE=/hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -syntenicNet \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 # Elapsed time: 400m22s cat fb.melGal1.chainTaeGut1Link.txt # 671105409 bases of 935922386 (71.705%) in intersection cd /hive/data/genomes/melGal1/bed ln -s lastzTaeGut1.2011-03-28 lastz.taeGut1 # running the swap mkdir /hive/data/genomes/taeGut1/bed/blastz.melGal1.swap cd /hive/data/genomes/taeGut1/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 36m7.017s cat fb.taeGut1.chainMelGal1Link.txt # 795001886 bases of 1222864691 (65.011%) in intersection cd /hive/data/genomes/taeGut1/bed ln -s blastz.melGal1.swap lastz.melGal1 ######################################################################### # SWAP hg19 lastz (DONE - 201i1-03-28 - Chin) # original alignment cat fb.hg19.chainMelGal1Link.txt # 76647912 bases of 2897316137 (2.645%) in intersection # and the swap mkdir /hive/data/genomes/melGal1/bed/blastz.hg19.swap cd /hive/data/genomes/melGal1/bed/blastz.hg19.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg19/bed/lastzMelGal1.2011-03-28/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 6m51.280s cat fb.melGal1.chainHg19Link.txt # 62120143 bases of 935922386 (6.637%) in intersection cd /hive/data/genomes/melGal1/bed ln -s blastz.hg19.swap lastz.hg19 ######################################################################### # SWAP mm9 lastz (DONE - 201i1-03-28 - Chin) # original alignment cat fb.mm9.chainMelGal1Link.txt # 62597891 bases of 2620346127 (2.389%) in intersection # running the swap mkdir /hive/data/genomes/melGal1/bed/blastz.mm9.swap cd /hive/data/genomes/melGal1/bed/blastz.mm9.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm9/bed/lastzMelGal1.2011-03-30/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 6m49.871s cat fb.melGal1.chainMm9Link.txt # 50649368 bases of 935922386 (5.412%) in intersection cd /hive/data/genomes/melGal1/bed ln -s blastz.mm9.swap lastz.mm9 ######################################################################### # SWAP monDom5 lastz (DONE - 2011-03-30 - Chin) # original alignment cat fb.monDom5.chainMelGal1Link.txt # 79077534 bases of 3501660299 (2.258%) in intersection # running the swap mkdir /hive/data/genomes/melGal1/bed/blastz.monDom5.swap cd /hive/data/genomes/melGal1/bed/blastz.monDom5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/monDom5/bed/lastzMelGal1.2011-03-30/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # Elapsed time: 10m15s cat fb.melGal1.chainMonDom5Link.txt # 64933013 bases of 935922386 (6.938%) in intersection cd /hive/data/genomes/melGal1/bed ln -s lastz.monDom5.swap lastz.monDom5 ######################################################################### # BLASTZ/CHAIN/NET Turkey vs Xenopus xenTro2 (DONE - 2011-04-02 - Chin) mkdir /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02 cd /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02 cat << '_EOF_' > DEF # Xenopus vs Turkey # TARGET: Turkey melGal1 SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit SEQ1_LEN=/scratch/data/melGal1/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 # QUERY: Xenopus xenTro2 - single chunk big enough to run entire chrom SEQ2_DIR=/scratch/data/xenTro2/xenTro2.2bit SEQ2_LEN=/scratch/data/xenTro2/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LAP=0 SEQ2_LIMIT=100 BASE=/hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 # real 288m3.884s cat fb.melGal1.chainXenTro2Link.txt # 36288270 bases of 935922386 (3.877%) in intersection cd /hive/data/genomes/melGal1/bed ln -s lastzXenTro2.2011-04-02 lastz.xenTro2 # running the swap mkdir /hive/data/genomes/xenTro2/bed/blastz.melGal1.swap cd /hive/data/genomes/xenTro2/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 3m40.138s cat fb.xenTro2.chainMelGal1Link.txt # 41967849 bases of 1359412157 (3.087%) in intersection cd /hive/data/genomes/xenTro2/bed ln -s blastz.melGal1.swap lastz.melGal1 ####################################################################### # lastz swap from danRer7 (DONE - 2011-04-02 - Chin) # original alignment to danRer7 cd /hive/data/genomes/danRer7/bed/lastzMelGal1.2011-04-01 cat fb.danRer7.chainMelGal1Link.txt # 35772051 bases of 1409770109 (2.537%) in intersection # running the swap mkdir /hive/data/genomes/melGal1/bed/blastz.danRer7.swap cd /hive/data/genomes/melGal1/bed/blastz.danRer7.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/danRer7/bed/lastzMelGal1.2011-04-01/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 3m28.672s cat fb.melGal1.chainDanRer7Link.txt # 25141503 bases of 935922386 (2.686%) in intersection cd /hive/data/genomes/melGal1/bed ln -s blastz.danRer7.swap lastz.danRer7 ####################################################################### # lastz swap from rn4 (DONE - 2011-04-04 - Chin) # original alignment to rn4 cd /hive/data/genomes/rn4/bed/lastzMelGal1.2011-04-04 cat fb.rn4.chainMelGal1Link.txt # 60383476 bases of 2571531505 (2.348%) in intersection # running the swap mkdir /hive/data/genomes/melGal1/bed/blastz.rn4.swap cd /hive/data/genomes/melGal1/bed/blastz.rn4.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn4/bed/lastzMelGal1.2011-04-04/DEF \ -swap \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 6m49.871s cat fb.melGal1.chainRn4Link.txt # 48636842 bases of 935922386 (5.197%) in intersection cd /hive/data/genomes/melGal1/bed ln -s blastz.rn4.swap lastz.rn4 ############################################################################ # HUMAN (hg18) PROTEINS TRACK (DONE 2011-04-26 braney) # bash if not using bash shell already cd /cluster/data/melGal1 mkdir /cluster/data/melGal1/blastDb awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst twoBitToFa -seqList=1meg.lst melGal1.2bit temp.fa faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft # 1140 pieces of 1140 written rm temp.fa 1meg.lst awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst twoBitToFa -seqList=less1meg.lst melGal1.2bit temp.fa faSplit about temp.fa 1000000 blastDb/y rm temp.fa less1meg.lst cd blastDb for i in *.fa do /hive/data/outside/blast229/formatdb -i $i -p F done rm *.fa ls *.nsq | wc -l # 1163 mkdir -p /cluster/data/melGal1/bed/tblastn.hg18KG cd /cluster/data/melGal1/bed/tblastn.hg18KG echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst wc -l query.lst # 1163 query.lst # we want around 250000 jobs calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(250000/`wc query.lst | awk '{print $1}'`\) # 36727/(250000/1163) = 170.854004 mkdir -p kgfa split -l 171 /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl kgfa/kg cd kgfa for i in *; do nice pslxToFa $i $i.fa; rm $i; done cd .. ls -1S kgfa/*.fa > kg.lst wc kg.lst # 215 215 2795 kg.lst mkdir -p blastOut for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done tcsh cd /cluster/data/melGal1/bed/tblastn.hg18KG cat << '_EOF_' > blastGsub #LOOP blastSome $(path1) {check in line $(path2)} {check out exists blastOut/$(root2)/q.$(root1).psl } #ENDLOOP '_EOF_' cat << '_EOF_' > blastSome #!/bin/sh BLASTMAT=/hive/data/outside/blast229/data export BLASTMAT g=`basename $2` f=/tmp/`basename $3`.$g for eVal in 0.01 0.001 0.0001 0.00001 0.000001 1E-09 1E-11 do if /hive/data/outside/blast229/blastall -M BLOSUM80 -m 0 -F no -e $eVal -p tblastn -d $1 -i $2 -o $f.8 then mv $f.8 $f.1 break; fi done if test -f $f.1 then if /cluster/bin/i386/blastToPsl $f.1 $f.2 then liftUp -nosort -type=".psl" -nohead $f.3 /cluster/data/melGal1/blastDb.lft carry $f.2 liftUp -nosort -type=".psl" -pslQ -nohead $3.tmp /cluster/data/hg18/bed/blat.hg18KG/protein.lft warn $f.3 if pslCheck -prot $3.tmp then mv $3.tmp $3 rm -f $f.1 $f.2 $f.3 $f.4 fi exit 0 fi fi rm -f $f.1 $f.2 $3.tmp $f.8 $f.3 $f.4 exit 1 '_EOF_' # << happy emacs chmod +x blastSome exit ssh swarm cd /cluster/data/melGal1/bed/tblastn.hg18KG gensub2 query.lst kg.lst blastGsub blastSpec para create blastSpec # para try, check, push, check etc. para time # Completed: 250045 of 250045 jobs # CPU time in finished jobs: 5922436s 98707.26m 1645.12h 68.55d 0.188 y # IO & Wait Time: 1789013s 29816.89m 496.95h 20.71d 0.057 y # Average job time: 31s 0.51m 0.01h 0.00d # Longest finished job: 330s 5.50m 0.09h 0.00d # Submission to last job: 12017s 200.28m 3.34h 0.14d ssh swarm cd /cluster/data/melGal1/bed/tblastn.hg18KG mkdir chainRun cd chainRun tcsh cat << '_EOF_' > chainGsub #LOOP chainOne $(path1) #ENDLOOP '_EOF_' cat << '_EOF_' > chainOne (cd $1; cat q.*.psl | simpleChain -prot -outPsl -maxGap=150000 stdin ../c.`basename $1`.psl) '_EOF_' chmod +x chainOne ls -1dS ../blastOut/kg?? > chain.lst gensub2 chain.lst single chainGsub chainSpec # do the cluster run for chaining para create chainSpec para try, check, push, check etc. # Completed: 215 of 215 jobs # CPU time in finished jobs: 1662s 27.70m 0.46h 0.02d 0.000 y # IO & Wait Time: 85567s 1426.12m 23.77h 0.99d 0.003 y # Average job time: 406s 6.76m 0.11h 0.00d # Longest finished job: 510s 8.50m 0.14h 0.01d # Submission to last job: 537s 8.95m 0.15h 0.01d cd /cluster/data/melGal1/bed/tblastn.hg18KG/blastOut for i in kg?? do cat c.$i.psl | awk "(\$13 - \$12)/\$11 > 0.6 {print}" > c60.$i.psl sort -rn c60.$i.psl | pslUniq stdin u.$i.psl awk "((\$1 / \$11) ) > 0.60 { print }" c60.$i.psl > m60.$i.psl echo $i done sort u.*.psl m60* | uniq | sort -T /tmp -k 14,14 -k 16,16n -k 17,17n > ../blastHg18KG.psl cd .. pslCheck blastHg18KG.psl # checked: 35263 failed: 0 errors: 0 # load table ssh hgwdev cd /cluster/data/melGal1/bed/tblastn.hg18KG hgLoadPsl melGal1 blastHg18KG.psl # check coverage featureBits melGal1 blastHg18KG # 18362443 bases of 935922386 (1.962%) in intersection featureBits melGal1 blastHg18KG ensGene -enrichment # blastHg18KG 1.962%, ensGene 2.680%, both 1.787%, cover 91.07%, enrich 33.99x rm -rf blastOut #end tblastn ############################################################################ # N-SCAN gene predictions (nscanGene) - (2011-07-07 markd) # obtained NSCAN predictions from Jeltje van Baren cd /hive/data/genomes/melGal1/bed/nscan cp -r /cluster/home/jeltje/hive/melGal1/for_browser/{melGal1.gtf, melGal1.prot.fa,readme.txt} . gzip melGal1.* chmod a-w * # load track gtfToGenePred -genePredExt melGal1.gtf.gz stdout | hgLoadGenePred -genePredExt melGal1 nscanGene stdin hgPepPred melGal1 generic nscanPep melGal1.prot.fa.gz rm *.tab # turkey/melGal1/trackDb.ra, add: track nscanGene override informant Nscan-EST on melGal1 with taeGut1 as informant and galGal3/taeGut1 parameter file. Not pseudogene masked. Used transmapped UCSC gene track as EST track. # verify top-level search spec, should produce no results and an error for egrep: hgsql -Ne 'select name from nscanGene' melGal1 | egrep -v '^chr[0-9a-zA-Z_]+\.([0-9]+|pasa)((\.[0-9a-z]+)?\.[0-9a-z]+)?$' |head ############################################################################ # lastz Turkey vs Frog xenTro3 (DONE - 2011-09-22 - Hiram) # This was repeated below with different lastz parameters that # give a bit higher coverage mkdir /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22 cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22 cat << '_EOF_' > DEF # Turkey vs. Frog # TARGET: Turkey melGal1 SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit SEQ1_LEN=/scratch/data/melGal1/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 # QUERY: Frog xenTro3 SEQ2_DIR=/scratch/data/xenTro3/xenTro3.2bit SEQ2_LEN=/scratch/data/xenTro3/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 SEQ2_LIMIT=30 BASE=/hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 & # real 648m6.498s cat fb.melGal1.chainXenTro3Link.txt # 35993896 bases of 935922386 (3.846%) in intersection cd /hive/data/genomes/melGal1/bed ln -s lastzXenTro3.2011-09-22 lastz.xenTro3 # running the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22/DEF \ -swap \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 3m29.407s cat fb.xenTro3.chainMelGal1Link.txt # 42384698 bases of 1358334882 (3.120%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.melGal1.swap lastz.melGal1 ####################################################################### # lastz Turkey vs Frog xenTro3 (DONE - 2011-09-27 - Hiram) mkdir /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27 cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27 cat << '_EOF_' > DEF # Turkey vs. Frog BLASTZ_H=2000 BLASTZ_Y=3400 BLASTZ_L=8000 BLASTZ_K=2200 BLASTZ_Q=/scratch/data/blastz/HoxD55.q # TARGET: Turkey melGal1 SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit SEQ1_LEN=/scratch/data/melGal1/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 # QUERY: Frog xenTro3 SEQ2_DIR=/scratch/data/xenTro3/xenTro3.2bit SEQ2_LEN=/scratch/data/xenTro3/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 SEQ2_LIMIT=30 BASE=/hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 & # real 537m35.283s cat fb.melGal1.chainXenTro3Link.txt # 49081729 bases of 935922386 (5.244%) in intersection cd /hive/data/genomes/melGal1/bed ln -s lastzXenTro3.2011-09-27 lastz.xenTro3 # running the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27/DEF \ -swap \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 4m9.023s cat fb.xenTro3.chainMelGal1Link.txt # 57087653 bases of 1358334882 (4.203%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.melGal1.swap lastz.melGal1 ####################################################################### # LASTZ Medium Ground Finch geoFor1 (DONE - 2012-07-28 - Hiram) # establish a screen to control this job with a name to indicate what it is screen -S melGal1 mkdir /hive/data/genomes/melGal1/bed/lastzGeoFor1.2012-07-28 cd /hive/data/genomes/melGal1/bed/lastzGeoFor1.2012-07-28 # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable # number of jobs, 50,000 to something under 100,000 cat << '_EOF_' > DEF # Turkey vs. Medium Ground Finch BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.02/bin/lastz # TARGET: turkey melGal1 SEQ1_DIR=/hive/data/genomes/melGal1/melGal1.2bit SEQ1_LEN=/hive/data/genomes/melGal1/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 SEQ1_LIMIT=70 # QUERY: Medium Ground Finch GeoFor1 SEQ2_DIR=/hive/data/genomes/geoFor1/geoFor1.2bit SEQ2_LEN=/hive/data/genomes/geoFor1/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 SEQ2_LIMIT=70 BASE=/hive/data/genomes/melGal1/bed/lastzGeoFor1.2012-07-28 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time doBlastzChainNet.pl -verbose=2 `pwd`/DEF \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ > do.log 2>&1 # Elapsed time: 917m34s cat fb.melGal1.chainGeoFor1Link.txt # 679715264 bases of 935922386 (72.625%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/melGal1/bed ln -s lastzGeoFor1.2012-07-28 lastz.geoFor1 mkdir /hive/data/genomes/geoFor1/bed/blastz.melGal1.swap cd /hive/data/genomes/geoFor1/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzGeoFor1.2012-07-28/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ > swap.log 2>&1 & # real 43m59.718s cat fb.geoFor1.chainHg19Link.txt # 698483400 bases of 1041286029 (67.079%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/geoFor1/bed ln -s blastz.melGal1.swap lastz.melGal1 ############################################################################## # lastz swap rn5 Rat (DONE - 2012-11-14 - Hiram) # the original alignment: cd /hive/data/genomes/rn5/bed/lastzMelGal1.2012-11-13 cat fb.rn5.chainMelGal1Link.txt # 93392251 bases of 2572853723 (3.630%) in intersection # and this swap: mkdir /hive/data/genomes/melGal1/bed/blastz.rn5.swap cd /hive/data/genomes/melGal1/bed/blastz.rn5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn5/bed/lastzMelGal1.2012-11-13/DEF \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ -swap -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 10m0.139s cat fb.melGal1.chainRn5Link.txt # 74638274 bases of 935922386 (7.975%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/melGal1/bed ln -s blastz.rn5.swap lastz.rn5 #########################################################################