src/hg/makeDb/doc/calJac3.txt 1.4
1.4 2010/02/18 22:15:04 hiram
trying to lastz papHam1 ogoGar1 micMur1 and tarSyr1 with difficulty
Index: src/hg/makeDb/doc/calJac3.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/calJac3.txt,v
retrieving revision 1.3
retrieving revision 1.4
diff -b -B -U 1000000 -r1.3 -r1.4
--- src/hg/makeDb/doc/calJac3.txt 16 Feb 2010 04:54:18 -0000 1.3
+++ src/hg/makeDb/doc/calJac3.txt 18 Feb 2010 22:15:04 -0000 1.4
@@ -1,606 +1,871 @@
# for emacs: -*- mode: sh; -*-
# $Id$
# Marmoset sequence: http://panda.genomics.org.cn/page/panda/download.jsp
# ftp.ncbi.nlm.nih.gov:genbank/genomes/Eukaryotes/vertebrates_mammals/
# Callithrix_jacchus/Callithrix_jacchus-3.2
# Callithrix jacchus
##########################################################################
# Download sequence (DONE - 2010-02-04 - Hiram)
mkdir /hive/data/genomes/calJac3
cd /hive/data/genomes/calJac3
mkdir genbank
cd genbank
wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \
--no-remove-listing -np \
"ftp.ncbi.nlm.nih.gov:genbank/genomes/Eukaryotes/vertebrates_mammals/Callithrix_jacchus/Callithrix_jacchus-3.2/*"
mkdir ucscChr
cd ucscChr
# fixup the accession names to become UCSC chrom names
zcat ../Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz \
| sed -e "s/^>.*gb|\([A-Z]*[0-9]*\).1.*/>chrUn_\1/" > chrUn.fa
zcat ../Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz \
| sed -e "s/^\([A-Z]*[0-9]*\).1/chrUn_\1/" > chrUn.agp
find ../Primary_Assembly/unlocalized_scaffolds/FASTA -type f \
| grep "unlocalized.scaf.fa.gz$" | head -1 | while read F
do
C=`basename ${F}`
C=${C/.unlocalized.scaf.fa.gz}
zcat "${F}" | sed -e "s/^>.*gb|\([A-Z]*[0-9]*\).1.*/>${C}_\1_random/"
done > chr_randoms.fa
find ../Primary_Assembly/unlocalized_scaffolds/AGP -type f | grep ".agp.gz$" \
| while read F
do
C=`basename ${F}`
echo -n "${C} "
C=${C/.unlocalized.scaf.agp.gz}
echo "${C}"
zcat "${F}" | sed -e "s/^\([A-Z]*[0-9]*\).1/${C}_\1_random/"
done > chr_randoms.agp
AC="../Primary_Assembly/assembled_chromosomes"
for F in ${AC}/FASTA/chr*.fa.gz
do
C=`basename $F`
C=${C/.fa.gz}
echo -n "${C} "
H=`zcat "${F}" | head -1`
chrN=`echo $H | sed -e "s/.*Callithrix jacchus chromosome /chr/; s/,
.*//"`
A=`echo $H | sed -e "s/. Callithrix.*//; s/.*gb.//"`
echo $chrN $A
zcat ${AC}/AGP/${chrN}.comp.agp.gz \
| sed -e "s/^${A}/${chrN}/" > ${chrN}.agp
echo ">${chrN}" > ${chrN}.fa
zcat ${AC}/FASTA/${chrN}.fa.gz | grep -v "^>" >> ${chrN}.fa
done
##########################################################################
# Initial genome build (DONE - 2009-12-17 - Hiram)
cd /hive/data/genomes/calJac3
cat << '_EOF_' > calJac3.config.ra
# Config parameters for makeGenomeDb.pl:
db calJac3
clade mammal
genomeCladePriority 16
scientificName Callithrix jacchus
commonName Marmoset
assemblyDate Mar. 2009
assemblyLabel WUSTL 3.2 (GCA_000004665.1)
orderKey 40
mitoAcc none
fastaFiles /hive/data/genomes/calJac3/genbank/ucscChr/*.fa
agpFiles /hive/data/genomes/calJac3/genbank/ucscChr/*.agp
# qualFiles none
dbDbSpeciesDir marmoset
taxId 9483
'_EOF_'
makeGenomeDb.pl -workhorse=hgwdev -stop=seq calJac3.config.ra > seq.out 2>&1
# real 4m5.924s
makeGenomeDb.pl -continue=agp -stop=agp calJac3.config.ra > agp.out 2>&1
# real 0m20.968s
makeGenomeDb.pl -continue=db -stop=db calJac3.config.ra > db.out 2>&1
# real 5m39.181s
XXX - chromInfo doesn't have large enough fields for the name keys
makeGenomeDb.pl -continue=dbDb -stop=dbDb calJac3.config.ra > dbDb.out 2>&1
makeGenomeDb.pl -continue=trackDb -stop=trackDb calJac3.config.ra > trackDb.out 2>&1
##########################################################################
# running repeat masker (DONE - 2010-02-02 - Hiram)
mkdir /hive/data/genomes/calJac3/bed/repeatMasker
cd /hive/data/genomes/calJac3/bed/repeatMasker
time doRepeatMasker.pl -buildDir=`pwd` -noSplit \
-bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \
-smallClusterHub=memk calJac3 > do.log 2>&1 &
# real 443m55.891s
cat faSize.rmsk.txt
# 2914958544 bases (162452744 N's 2752505800 real 1439949519 upper 1312556281
# lower) in 14205 sequences in 1 files
# %45.03 masked total, %47.69 masked real
##########################################################################
# running simple repeat (DONE - 2010-02-02 - Hiram)
mkdir /hive/data/genomes/calJac3/bed/simpleRepeat
cd /hive/data/genomes/calJac3/bed/simpleRepeat
time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \
-dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \
calJac3 > do.log 2>&1 &
# real 198m33.953s
cat fb.simpleRepeat
# 66211508 bases of 2752505800 (2.405%) in intersection
cd /hive/data/genomes/calJac3
twoBitMask calJac3.rmsk.2bit \
-add bed/simpleRepeat/trfMask.bed calJac3.2bit
# you can safely ignore the warning about fields >= 13
twoBitToFa calJac3.2bit stdout | faSize stdin > faSize.calJac3.2bit.txt
cat faSize.calJac3.2bit.txt
# 2914958544 bases (162452744 N's 2752505800 real 1439244378 upper
# 1313261422 lower) in 14205 sequences in 1 files
# %45.05 masked total, %47.71 masked real
rm /gbdb/calJac3/calJac3.2bit
ln -s `pwd`/calJac3.2bit /gbdb/calJac3/calJac3.2bit
#########################################################################
# MAKE 11.OOC FILES FOR BLAT (DONE - 2010-02-03 - Hiram)
ssh kolossus
# numerator is calJac3 gapless bases "real" as reported by faSize
# denominator is hg17 gapless bases as reported by featureBits,
# 1024 is threshold used for human -repMatch:
calc \( 2752505800 / 2897310462 \) \* 1024
# ( 2752505800 / 2897310462 ) * 1024 = 972.821510
# ==> use -repMatch=950 according to size scaled down from 1024 for human.
# and rounded down to nearest 50
cd /hive/data/genomes/calJac3
blat calJac3.2bit /dev/null /dev/null -tileSize=11 \
-makeOoc=jkStuff/calJac3.11.ooc -repMatch=950
# Wrote 32908 overused 11-mers to jkStuff/calJac3.11.ooc
mkdir /hive/data/staging/data/calJac3
cp -p calJac3.2bit chrom.sizes jkStuff/calJac3.11.ooc \
/hive/data/staging/data/calJac3
gapToLift -bedFile=jkStuff/nonBridgedGaps.bed calJac3 \
jkStuff/calJac3.nonBridged.lft
##########################################################################
# BLATSERVERS ENTRY (DONE - 2009-12-23 - Hiram)
# After getting a blat server assigned by the Blat Server Gods,
ssh hgwdev
hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
VALUES ("calJac3", "blat12", "17804", "1", "0"); \
INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
VALUES ("calJac3", "blat12", "17805", "0", "1");' \
hgcentraltest
# test it with some sequence
############################################################################
# reset position to RHO location as found from blat of hg19 RHO gene
hgsql -e \
'update dbDb set defaultPos="chr15:67439781-67471107" where name="calJac3";' \
hgcentraltest
############################################################################
# genbank run (DONE - 2010-02-04 - Hiram)
ssh hgwdev
cd $HOME/kent/src/hg/makeDb/genbank
# edit etc/genbank.conf to add this section just before calJac1:
# Marmoset
calJac3.serverGenome = /hive/data/genomes/calJac3/calJac3.2bit
calJac3.clusterGenome = /scratch/data/calJac3/calJac3.2bit
calJac3.ooc = /scratch/data/calJac3/calJac3.11.ooc
calJac3.lift = no
calJac3.perChromTables = no
calJac3.refseq.mrna.native.pslCDnaFilter = ${ordered.refseq.mrna.native.pslCDnaFilter}
calJac3.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter}
calJac3.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter}
calJac3.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter}
calJac3.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter}
calJac3.genbank.est.xeno.pslCDnaFilter = ${ordered.genbank.est.xeno.pslCDnaFilter}
calJac3.downloadDir = calJac3
calJac3.refseq.mrna.native.load = yes
calJac3.refseq.mrna.xeno.load = yes
calJac3.refseq.mrna.xeno.loadDesc = yes
cvs ci -m "adding marmoset calJac3" etc/genbank.conf
make etc-update
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 calJac3 &
# var/build/logs/2010.02.10-09:19:20.calJac3.initalign.log
# real 238m59.390s
ssh hgwdev
cd /cluster/data/genbank
time ./bin/gbDbLoadStep -drop -initialLoad calJac3 &
# logFile: var/dbload/hgwdev/logs/2010.02.11-10:41:22.dbload.log
# real 25m55.474s
# enable daily alignment and update of hgwdev
cd ~/kent/src/hg/makeDb/genbank
cvsup
# add calJac3 to:
etc/align.dbs
etc/hgwdev.dbs
cvs ci -m "Adding calJac3 - Marmoset - Callithrix jacchus" \
etc/align.dbs etc/hgwdev.dbs
make etc-update
# done - 2010-02-11 - Hiram
############################################################################
# running cpgIsland business (DONE - 2010-02-11 - Hiram)
mkdir /hive/data/genomes/calJac3/bed/cpgIsland
cd /hive/data/genomes/calJac3/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 ../../calJac3.2bit:$C stdout \
| maskOutFa stdin hard hardMaskedFa/${C}.fa
done
ssh swarm
cd /hive/data/genomes/calJac3/bed/cpgIsland
mkdir results
cut -f1 ../../chrom.sizes > chr.list
cat << '_EOF_' > template
#LOOP
./runOne $(root1) {check 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
gensub2 chr.list single template jobList
para create jobList
para try
para check ... etc
para time
# Completed: 14205 of 14205 jobs
# CPU time in finished jobs: 205s 3.41m 0.06h 0.00d 0.000 y
# IO & Wait Time: 38701s 645.02m 10.75h 0.45d 0.001 y
# Average job time: 3s 0.05m 0.00h 0.00d
# Longest finished job: 23s 0.38m 0.01h 0.00d
# Submission to last job: 202s 3.37m 0.06h 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
cd /hive/data/genomes/calJac3/bed/cpgIsland
hgLoadBed calJac3 cpgIslandExt -tab \
-sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
# Reading cpgIsland.bed
# Loaded 32732 elements of size 10
# Sorted
# Saving bed.tab
# Loading calJac3
# cleanup
rm -fr hardMaskedFa
############################################################################
# LASTZ Human Swap (DONE - 2010-02-11 - Hiram)
# original alignment on hg19:
cd /hive/data/genomes/hg19/bed/lastzCalJac3.2010-02-11
cat fb.hg19.chainCalJac3Link.txt
# 2047068864 bases of 2897316137 (70.654%) in intersection
# and for this swap
mkdir /hive/data/genomes/calJac3/bed/blastz.hg19.swap
cd /hive/data/genomes/calJac3/bed/blastz.hg19.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/hg19/bed/lastzCalJac3.2010-02-11/DEF \
-swap -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
# real 120m42.991s
cat fb.calJac3.chainHg19Link.txt
# 2030475813 bases of 2752505800 (73.768%) in intersection
#####################################################################
# LASTZ Rhesus Swap (DONE - 2010-02-11 - Hiram)
# original alignment to rheMac2
cd /hive/data/genomes/rheMac2/bed/lastzCalJac3.2010-02-11
cat fb.rheMac2.chainCalJac3Link.txt
# 1871513554 bases of 2646704109 (70.711%) in intersection
# and for this swap
mkdir /hive/data/genomes/calJac3/bed/blastz.rheMac2.swap
cd /hive/data/genomes/calJac3/bed/blastz.rheMac2.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/rheMac2/bed/lastzCalJac3.2010-02-11/DEF \
-swap -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
# real 142m34.894s
cat fb.calJac3.chainHg19Link.txt
# 1916431926 bases of 2752505800 (69.625%) in intersection
############################################################################
# LASTZ Chimp Swap (DONE - 2010-02-11 - Hiram)
# original alignment to panTro2
cd /hive/data/genomes/panTro2/bed/lastzCalJac3.2010-02-11
cat fb.panTro2.chainCalJac3Link.txt
# 2016331285 bases of 2909485072 (69.302%) in intersection
# and this swap run
mkdir /hive/data/genomes/calJac3/bed/blastz.panTro2.swap
cd /hive/data/genomes/calJac3/bed/blastz.panTro2.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/panTro2/bed/lastzCalJac3.2010-02-11/DEF \
-swap -noLoadChainSplit -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
# real 118m42.203s
cat fb.calJac3.chainHg19Link.txt
# 1990168262 bases of 2752505800 (72.304%) in intersection
############################################################################
# LASTZ Orangutan Swap (DONE - 2010-02-11 - Hiram)
# original alignment to ponAbe2
cd /hive/data/genomes/ponAbe2/bed/lastzCalJac3.2010-02-11
cat fb.ponAbe2.chainCalJac3Link.txt
# 2086557592 bases of 3093572278 (67.448%) in intersection
# and this swap run
mkdir /hive/data/genomes/calJac3/bed/blastz.ponAbe2.swap
cd /hive/data/genomes/calJac3/bed/blastz.ponAbe2.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/ponAbe2/bed/lastzCalJac3.2010-02-11/DEF \
-swap -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
# real 146m12.301s
cat fb.calJac3.chainHg19Link.txt
# 1978857628 bases of 2752505800 (71.893%) in intersection
#####################################################################
# lastz Dog CanFam2 Swap (DONE - 2010-02-14 - Hiram)
# original alignment
cd /hive/data/genomes/canFam2/bed/lastzCalJac3.2010-02-12
cat fb.canFam2.chainCalJac3Link.txt
# 1363307334 bases of 2384996543 (57.162%) in intersection
# and for the swap
mkdir /hive/data/genomes/calJac3/bed/blastz.canFam2.swap
cd /hive/data/genomes/calJac3/bed/blastz.canFam2.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/canFam2/bed/lastzCalJac3.2010-02-12/DEF \
-swap -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 129m56.144s
cat fb.calJac3.chainHg19Link.txt
# 1397333116 bases of 2752505800 (50.766%) in intersection
#########################################################################
# lastz Mouse Mm9 Swap (DONE - 2010-02-14 - Hiram)
# original alignment to mouse
cd /hive/data/genomes/mm9/bed/lastzCalJac3.2010-02-12
cat fb.mm9.chainCalJac3Link.txt
# 859869647 bases of 2620346127 (32.815%) in intersection
# and for the swap
mkdir /hive/data/genomes/calJac3/bed/blastz.mm9.swap
cd /hive/data/genomes/calJac3/bed/blastz.mm9.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/mm9/bed/lastzCalJac3.2010-02-12/DEF \
-swap -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 90m38.739s
cat fb.calJac3.chainHg19Link.txt
# 861811978 bases of 2752505800 (31.310%) in intersection
#####################################################################
# lastz Opossum monDom5 Swap (DONE - 2010-02-14 - Hiram)
# original alignment to Opossum
cd /hive/data/genomes/monDom5/bed/lastzCalJac3.2010-02-12
cat fb.monDom5.chainCalJac3Link.txt
# 216197506 bases of 3501660299 (6.174%) in intersection
# and for the swap
mkdir /hive/data/genomes/calJac3/bed/blastz.monDom5.swap
cd /hive/data/genomes/calJac3/bed/blastz.monDom5.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/monDom5/bed/lastzCalJac3.2010-02-12/DEF \
-swap -noLoadChainSplit -syntenicNet \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 &
# real 110m13.435s
cat fb.calJac3.chainMonDom5Link.txt
# 217614612 bases of 2752505800 (7.906%) in intersection
##############################################################################
# HUMAN (hg18) PROTEINS TRACK (working braney...)
# bash if not using bash shell already
cd /cluster/data/calJac3
mkdir /cluster/data/calJac3/blastDb
awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
twoBitToFa -seqList=1meg.lst calJac3.2bit temp.fa
faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
rm temp.fa 1meg.lst
awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst
twoBitToFa -seqList=less1meg.lst calJac3.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
# 3275
mkdir -p /cluster/data/calJac3/bed/tblastn.hg18KG
cd /cluster/data/calJac3/bed/tblastn.hg18KG
echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst
wc -l query.lst
# 3275 query.lst
# we want around 350000 jobs
calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(350000/`wc query.lst | awk '{print $1}'`\)
# 36727/(350000/3275) = 343.659786
mkdir -p kgfa
split -l 344 /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
# 107 107 1391 kg.lst
mkdir -p blastOut
for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done
tcsh
cd /cluster/data/calJac3/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/calJac3/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/calJac3/bed/tblastn.hg18KG
gensub2 query.lst kg.lst blastGsub blastSpec
para create blastSpec
# para try, check, push, check etc.
para time
# Completed: 100076 of 100076 jobs
# CPU time in finished jobs: 2219419s 36990.31m 616.51h 25.69d 0.070 y
# IO & Wait Time: 445463s 7424.39m 123.74h 5.16d 0.014 y
# Average job time: 27s 0.44m 0.01h 0.00d
# Longest finished job: 319s 5.32m 0.09h 0.00d
# Submission to last job: 2791s 46.52m 0.78h 0.03d
ssh swarm
cd /cluster/data/calJac3/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=12000 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: 254 of 254 jobs
# CPU time in finished jobs: 590077s 9834.62m 163.91h 6.83d 0.019 y
# IO & Wait Time: 20749s 345.81m 5.76h 0.24d 0.001 y
# Average job time: 2405s 40.08m 0.67h 0.03d
# Longest finished job: 36592s 609.87m 10.16h 0.42d
# Submission to last job: 36604s 610.07m 10.17h 0.42d
cd /cluster/data/calJac3/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: 38517 failed: 0 errors: 0
# load table
ssh hgwdev
cd /cluster/data/calJac3/bed/tblastn.hg18KG
hgLoadPsl calJac3 blastHg18KG.psl
# check coverage
featureBits calJac3 blastHg18KG
# 18646943 bases of 332311746 (5.611%) in intersection
featureBits calJac3 blastHg18KG ensGene -enrichment
# blastHg18KG 5.611%, ensGene 9.520%, both 4.993%, cover 88.97%, enrich 9.35x
rm -rf blastOut
#end tblastn
##############################################################################
+# papHam1 Baboon LASTZ/CHAIN/NET (DONE - 2010-02-15 - Hiram)
+ screen # use a screen to manage this multi-day job
+ mkdir /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15
+ cd /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15
+
+ cat << '_EOF_' > DEF
+# baboon vs. marmoset
+# same paramters as human hg19 vs marmoset calJac3
+BLASTZ=lastz
+# maximum M allowed with lastz is only 254
+BLASTZ_M=254
+BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q
+BLASTZ_O=600
+BLASTZ_E=150
+BLASTZ_K=4500
+BLASTZ_Y=15000
+BLASTZ_T=2
+
+# TARGET: Marmoset (calJac3)
+SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit
+SEQ1_LEN=/scratch/data/calJac3/chrom.sizes
+SEQ1_LIMIT=50
+SEQ1_CHUNK=20000000
+SEQ1_LAP=10000
+
+# QUERY: Baboon papHam1
+SEQ2_DIR=/scratch/data/papHam1/papHam1.2bit
+SEQ2_LEN=/scratch/data/papHam1/chrom.sizes
+SEQ2_CHUNK=20000000
+SEQ2_LIMIT=300
+SEQ2_LAP=0
+
+BASE=/hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15
+TMPDIR=/scratch/tmp
+'_EOF_'
+ # << this line keeps emacs coloring happy
+
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+XXX - running Mon Feb 15 21:49:52 PST 2010
+ # real 287m24.258s
+ cat fb.calJac3.chainPapHam1Link.txt
+ # 2047068864 bases of 2897316137 (70.654%) in intersection
+
+ mkdir /hive/data/genomes/papHam1/bed/blastz.calJac3.swap
+ cd /hive/data/genomes/papHam1/bed/blastz.calJac3.swap
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-11/DEF \
+ -swap -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
+ # real 120m42.991s
+ cat fb.papHam1.chainCalJac3Link.txt
+ # 2030475813 bases of 2752505800 (73.768%) in intersection
+
+##############################################################################
+# tarSyr1 Tarsier LASTZ/CHAIN/NET (DONE - 2010-02-15 - Hiram)
+ screen # use a screen to manage this multi-day job
+ mkdir /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-15
+ cd /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-15
+
+ cat << '_EOF_' > DEF
+# tarsier vs. marmoset
+# same paramters as human hg19 vs marmoset calJac3
+BLASTZ=lastz
+# maximum M allowed with lastz is only 254
+BLASTZ_M=254
+BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q
+BLASTZ_O=600
+BLASTZ_E=150
+BLASTZ_K=4500
+BLASTZ_Y=15000
+BLASTZ_T=2
+
+# TARGET: Marmoset (calJac3)
+SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit
+SEQ1_LEN=/scratch/data/calJac3/chrom.sizes
+SEQ1_LIMIT=50
+SEQ1_CHUNK=20000000
+SEQ1_LAP=10000
+
+# QUERY: Tarsier tarSyr1
+SEQ2_DIR=/scratch/data/tarSyr1/tarSyr1.2bit
+SEQ2_LEN=/scratch/data/tarSyr1/chrom.sizes
+SEQ2_CHUNK=20000000
+SEQ2_LIMIT=300
+SEQ2_LAP=0
+
+BASE=/hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-15
+TMPDIR=/scratch/tmp
+'_EOF_'
+ # << this line keeps emacs coloring happy
+
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+XXX - running Mon Feb 15 22:00:23 PST 2010
+ # real 287m24.258s
+ cat fb.calJac3.chainTarSyr1Link.txt
+ # 2047068864 bases of 2897316137 (70.654%) in intersection
+
+ mkdir /hive/data/genomes/tarSyr1/bed/blastz.calJac3.swap
+ cd /hive/data/genomes/tarSyr1/bed/blastz.calJac3.swap
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-11/DEF \
+ -swap -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
+ # real 120m42.991s
+ cat fb.tarSyr1.chainCalJac3Link.txt
+ # 2030475813 bases of 2752505800 (73.768%) in intersection
+
+#####################################################################
+# micMur1 Mouse lemur LASTZ/CHAIN/NET (DONE - 2010-02-15 - Hiram)
+# Mouse lemur ( Microcebus murinus)
+ screen # use a screen to manage this multi-day job
+ mkdir /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-15
+ cd /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-15
+
+ cat << '_EOF_' > DEF
+# mouse lemur vs. marmoset
+# same paramters as human hg19 vs marmoset calJac3
+BLASTZ=lastz
+# maximum M allowed with lastz is only 254
+BLASTZ_M=254
+BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q
+BLASTZ_O=600
+BLASTZ_E=150
+BLASTZ_K=4500
+BLASTZ_Y=15000
+BLASTZ_T=2
+
+# TARGET: Marmoset (calJac3)
+SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit
+SEQ1_LEN=/scratch/data/calJac3/chrom.sizes
+SEQ1_LIMIT=50
+SEQ1_CHUNK=20000000
+SEQ1_LAP=10000
+
+# QUERY: Mouse lemur micMur1
+SEQ2_DIR=/hive/data/genomes/micMur1/micMur1.2bit
+SEQ2_LEN=/hive/data/genomes/micMur1/chrom.sizes
+SEQ2_CHUNK=20000000
+SEQ2_LIMIT=300
+SEQ2_LAP=0
+
+BASE=/hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-15
+TMPDIR=/scratch/tmp
+'_EOF_'
+ # << this line keeps emacs coloring happy
+
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+BASE=/hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-17
+XXX - running Wed Feb 17 13:57:47 PST 2010
+
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+ # failed first kluster job, finished manually
+ # real 287m24.258s
+ cat fb.calJac3.chainMicMur1Link.txt
+ # 2047068864 bases of 2897316137 (70.654%) in intersection
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ -continue=cat `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > cat.log 2>&1 &
+XXX - running Wed Feb 17 10:36:56 PST 2010
+
+ mkdir /hive/data/genomes/micMur1/bed/blastz.calJac3.swap
+ cd /hive/data/genomes/micMur1/bed/blastz.calJac3.swap
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-11/DEF \
+ -swap -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
+ # real 120m42.991s
+ cat fb.micMur1.chainCalJac3Link.txt
+ # 2030475813 bases of 2752505800 (73.768%) in intersection
+
+#####################################################################
+# otoGar1 Bushbaby LASTZ/CHAIN/NET (DONE - 2010-02-15 - Hiram)
+ screen # use a screen to manage this multi-day job
+ mkdir /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-15
+ cd /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-15
+
+ cat << '_EOF_' > DEF
+# bushbaby vs. marmoset
+# same paramters as human hg19 vs marmoset calJac3
+BLASTZ=lastz
+# maximum M allowed with lastz is only 254
+BLASTZ_M=254
+BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q
+BLASTZ_O=600
+BLASTZ_E=150
+BLASTZ_K=4500
+BLASTZ_Y=15000
+BLASTZ_T=2
+
+# TARGET: Marmoset (calJac3)
+SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit
+SEQ1_LEN=/scratch/data/calJac3/chrom.sizes
+SEQ1_LIMIT=50
+SEQ1_CHUNK=20000000
+SEQ1_LAP=10000
+
+# QUERY: Bushbaby otoGar1
+SEQ2_DIR=/scratch/data/otoGar1/otoGar1.rmsk.2bit
+SEQ2_LEN=/scratch/data/otoGar1/chrom.sizes
+SEQ2_CHUNK=20000000
+SEQ2_LIMIT=300
+SEQ2_LAP=0
+
+BASE=/hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-15
+TMPDIR=/scratch/tmp
+'_EOF_'
+ # << this line keeps emacs coloring happy
+
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+BASE=/hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-17
+XXX - running Wed Feb 17 13:57:47 PST 2010
+
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > do.log 2>&1 &
+ # failed first kluster job, finished manually
+ # real 287m24.258s
+ cat fb.calJac3.chainOtoGar1Link.txt
+ # 2047068864 bases of 2897316137 (70.654%) in intersection
+ time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \
+ -continue=cat `pwd`/DEF \
+ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ > cat.log 2>&1 &
+XXX - running Wed Feb 17 10:36:21 PST 2010
+
+ mkdir /hive/data/genomes/otoGar1/bed/blastz.calJac3.swap
+ cd /hive/data/genomes/otoGar1/bed/blastz.calJac3.swap
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-11/DEF \
+ -swap -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
+ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 &
+ # real 120m42.991s
+ cat fb.otoGar1.chainCalJac3Link.txt
+ # 2030475813 bases of 2752505800 (73.768%) in intersection
+
+#####################################################################