src/hg/makeDb/doc/bosTau4.txt 1.19
1.19 2009/11/25 21:48:37 hiram
change autoScaleDefault to autoScale
Index: src/hg/makeDb/doc/bosTau4.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/bosTau4.txt,v
retrieving revision 1.18
retrieving revision 1.19
diff -b -B -U 1000000 -r1.18 -r1.19
--- src/hg/makeDb/doc/bosTau4.txt 20 Sep 2009 17:16:41 -0000 1.18
+++ src/hg/makeDb/doc/bosTau4.txt 25 Nov 2009 21:48:37 -0000 1.19
@@ -1,1816 +1,1816 @@
# for emacs: -*- mode: sh; -*-
# Bos Taurus -- Baylor Release 4.0 (Oct 4 2007)
# "$Id$"
#########################################################################
# DOWNLOAD SEQUENCE (DONE - 2008-03-07 - Hiram)
ssh kkstore06
mkdir /cluster/store3/bosTau4
ln -s /cluster/store3/bosTau4 /cluster/data/bosTau4
mkdir /cluster/data/bosTau4/baylor
cd /cluster/data/bosTau4/baylor
for F in README.Btau20070913.txt Contigs/*
do
wget --timestamping \
"ftp://ftp.hgsc.bcm.tmc.edu/pub/data/Btaurus/fasta/Btau20070913-freeze/${F}"
done
mkdir chroms
cd chroms
wget --timestamping \
"ftp://ftp.hgsc.bcm.tmc.edu/pub/data/Btaurus/fasta/Btau20070913-freeze/LinearScaffolds/*"
# Looking up chrM in Entrez nucleotide search,
# from bosTau3 note: NC_006853 GI:60101824
# search for Bos taurus mitochondrion complete genome
# Note, their Contigs files are from their 2006 release.
# Their chrom files are new to this release.
# fixup the Chr case names in the agp file:
sed -e "s/^Chr/chr/" Btau20070913.All.ac.agp > bosTau4.agp
# fixup the contig names in the contigs.fa and qual files, utilize
# the following perl scripts:
cat << '_EOF_' > fixContigNames.pl
#!/usr/bin/env perl
use warnings;
use strict;
my $argc = scalar(@ARGV);
if (1 != $argc) {
print "usage ./fixContigNames.pl Btau20060815.contigs.fa.gz \\\n",
" | gzip > bosTau4.contigs.fa.gz\n";
exit 255;
}
my $fName = shift;
open (FH,"zcat $fName|") or die "can not read $fName $!";
while (my $line = <FH>) {
if ($line =~ m/^>/) {
my @a = split('\|', $line);
$a[3] =~ s/\.[0-9]+//;
print ">$a[3]\n";
} else {
print "$line";
}
}
close (FH)
'_EOF_'
# << happy emacs
cat << '_EOF_' > fixQuals.pl
#!/usr/bin/env perl
use warnings;
use strict;
my $argc = scalar(@ARGV);
if (1 != $argc) {
print "usage ./fixQuals.pl Btau20060815.contigs.fa.qual.gz \\\n",
" | gzip > bosTau4.qual.gz\n";
exit 255;
}
my $fName = shift;
open (FH,"zcat $fName|") or die "can not read $fName $!";
while (my $line = <FH>) {
if ($line =~ m/^>/) {
$line =~ s/\.[0-9]+.*//;
}
print "$line";
}
close (FH)
'_EOF_'
# << happy emacs
# There are extra records in the qual and contigs fa files.
# Only the sequences in the agp file should be used:
grep "^chr" bosTau4.agp | egrep -v "clone|fragment" | cut -f6 \
| sort > agp.names
# run those scripts, extract only the used records, and lift the
# qual scores via the AGP file to a qac file.
chmod +x fixContigNames.pl fixQuals.pl
./fixContigNames.pl Btau20060815.contigs.fa.gz \
| faSomeRecords stdin agp.names stdout | gzip > bosTau4.contigs.fa.gz
./fixQuals.pl Btau20060815.contigs.fa.qual.gz \
| faSomeRecords stdin agp.names stdout \
| qaToQac stdin stdout \
| qacAgpLift bosTau4.agp stdin bosTau4.qual.qac
#############################################################################
# running makeGenomeDb.pl (DONE - 2008-03-07 - Hiram)
ssh kkstore06
cd /cluster/data/bosTau4
cat << '_EOF_' > bosTau4.config.ra
db bosTau4
clade mammal
scientificName Bos Taurus
assemblyDate Oct. 2007
assemblyLabel Baylor Release 4.0
orderKey 236
dbDbSpeciesDir cow
mitoAcc 60101824
agpFiles /cluster/data/bosTau4/baylor/bosTau4.agp
fastaFiles /cluster/data/bosTau4/baylor/bosTau4.contigs.fa.gz
qualFiles /cluster/data/bosTau4/baylor/bosTau4.qual.qac
commonName Cow
'_EOF_'
# << happy emacs
makeGenomeDb.pl bosTau4.config.ra > makeGenomeDb.log 2>&1
#########################################################################
# REPEATMASKER (DONE - 2008-03-07 - Hiram)
ssh kkstore06
screen # use a screen to manage this job
mkdir /cluster/data/bosTau4/bed/repeatMasker
cd /cluster/data/bosTau4/bed/repeatMasker
#
# This was going too slow on memk, on the order of 8 days.
# chilled that memk batch, then switched it to finish on pk
doRepeatMasker.pl -buildDir=/cluster/data/bosTau4/bed/repeatMasker \
-bigClusterHub=memk bosTau4 > do.log 2>&1 &
# continuing
doRepeatMasker.pl -buildDir=/cluster/data/bosTau4/bed/repeatMasker \
-continue=cat -bigClusterHub=pk bosTau4 > cat.log 2>&1 &
time nice -n +19 featureBits bosTau4 rmsk > fb.bosTau4.rmsk.txt 2>&1 &
# 1280927549 bases of 2731830700 (46.889%) in intersection
# RepeatMasker and lib version from do.log:
# RepeatMasker,v 1.20 2008/01/16 18:15:45 angie Exp $
# Jan 11 2008 (open-3-1-9) version of RepeatMasker
# CC RELEASE 20071204;
# Compare coverage to previous assembly:
featureBits bosTau3 rmsk
# 1200525422 bases of 2731807384 (43.946%) in intersection
featureBits bosTau2 rmsk
# 1291561904 bases of 2812203870 (45.927%) in intersection
#########################################################################
# SIMPLE REPEATS TRF (DONE - 2008-03-07 - Hiram)
ssh kkstore06
screen # use a screen to manage this job
mkdir /cluster/data/bosTau4/bed/simpleRepeat
cd /cluster/data/bosTau4/bed/simpleRepeat
#
doSimpleRepeat.pl -buildDir=/cluster/data/bosTau4/bed/simpleRepeat \
bosTau4 > do.log 2>&1 &
# after RM run is done, add this mask:
cd /cluster/data/bosTau4
twoBitMask bosTau4.rmsk.2bit -add bed/simpleRepeat/trfMask.bed bosTau4.2bit
twoBitToFa bosTau4.2bit stdout | faSize stdin
# 2917974530 bases (186161294 N's 2731813236 real
# 1450855409 upper 1280957827 lower) in 11900 sequences in 1 files
# %43.90 masked total, %46.89 masked real
twoBitToFa bosTau4.rmsk.2bit stdout | faSize stdin
# 2917974530 bases (186161294 N's 2731813236 real
# 1451369861 upper 1280443375 lower) in 11900 sequences in 1 files
# %43.88 masked total, %46.87 masked real
#########################################################################
# Create OOC file for genbank runs (DONE - 2008-03-10 - Hiram)
# use same repMatch value as bosTau2
ssh kkstore06
cd /cluster/data/bosTau3
blat bosTau4.2bit /dev/null /dev/null -tileSize=11 \
-makeOoc=bosTau4.11.ooc -repMatch=1005
ssh kkr1u00
mkdir /iscratch/i/bosTau4
cd /iscratch/i/bosTau4
cp -p /cluster/data/bosTau4/bosTau4.2bit .
for R in 2 3 4 5 6 7 8
do
rsync -a --progress ./ kkr${R}u00:/iscratch/i/bosTau4/
done
#########################################################################
# Starting Genbank
ssh hgwdev
cd $HOME/kent/src/hg/makeDb/genbank
# edit etc/genbank.conf to add the following entry:
# bosTau4 (B. taurus) 31 chromosomes plus 11869 scaffolds
bosTau4.serverGenome = /cluster/data/bosTau4/bosTau4.2bit
bosTau4.clusterGenome = /iscratch/i/bosTau4/bosTau4.2bit
bosTau4.ooc = /cluster/data/bosTau4/bosTau4.11.ooc
bosTau4.refseq.mrna.native.pslCDnaFilter = ${lowCover.refseq.mrna.native.pslCDnaFilter}
bosTau4.refseq.mrna.xeno.pslCDnaFilter = ${lowCover.refseq.mrna.xeno.pslCDnaFilter}
bosTau4.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter}
bosTau4.genbank.mrna.xeno.pslCDnaFilter = ${lowCover.genbank.mrna.xeno.pslCDnaFilter}
bosTau4.genbank.est.native.pslCDnaFilter = ${lowCover.genbank.est.native.pslCDnaFilter}
bosTau4.lift = no
bosTau4.downloadDir = bosTau4
bosTau4.perChromTables = no
bosTau4.mgc = yes
cvs ci -m "Added bosTau4." etc/genbank.conf
# update /cluster/data/genbank/:
make etc-update
ssh genbank
screen # control this business with a screen since it takes a while
cd /cluster/data/genbank
# This is a call to a script that will push our jobs out to the cluster
# since it's a big job.
time nice -n +19 bin/gbAlignStep -initial bosTau4 &
# logFile: var/build/logs/2008.03.10-14:14:43.bosTau4.initalign.log
# real 567m7.431s
# load database when finished
ssh hgwdev
cd /cluster/data/genbank
time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad bosTau4
# logFile: var/dbload/hgwdev/logs/2008.03.11-08:48:17.dbload.log
# real 34m55.565s
# enable daily alignment and update of hgwdev (DONE - 2008-03-11 - Hiram)
cd ~/kent/src/hg/makeDb/genbank
cvsup
# add bosTau4 to:
etc/align.dbs
etc/hgwdev.dbs
cvs ci -m "Added bosTau4." etc/align.dbs etc/hgwdev.dbs
make etc-update
############################################################################
# DONE - 2008-03-11 - Hiram
# Reset default position to an area of chr3 with a number of genes
hgsql -e \
'update dbDb set defaultPos="chr3:21083267-21232422" where name="bosTau4";' \
hgcentraltest
# And there was a mistake in the description date entry
hgsql -e \
'update dbDb set description="Oct. 2007" where name="bosTau4";' \
hgcentraltest
#########################################################################
## genscan run (DONE - 2008-03-11 - Hiram)
## create hard masked sequence
ssh kkstore06
cd /cluster/data/bosTau4
mkdir hardMasked
for C in `cut -f1 chrom.sizes`
do
echo "hardMasked/${C}.hard.fa"
twoBitToFa -seq=${C} bosTau4.2bit stdout \
| maskOutFa stdin hard hardMasked/${C}.hard.fa
ls -ld "hardMasked/${C}.hard.fa"
done
# And, make sure there aren't any sequences in this lot that have
# become all N's with no sequence left in them. This drives genscan nuts
echo hardMasked/*.hard.fa | xargs faCount > faCount.hard.txt
# the lowest three are:
egrep -v "^#|^total" faCount.hard.txt \
| awk '{print $1,$2-$7}' | sort -k2,2nr | tail -3
# chrUn.004.9957 0
# chrUn.004.9961 0
# chrUn.004.9996 0
# There are a whole bunch of these, and many with just a few bases.
# Actually, before removing these for genscan, run the cpgIsland
# business first since it can work on them all.
# So, remove any with less than 100 bases of sequence
egrep -v "^#|^total" faCount.hard.txt \
| awk '{size=$2-$7; if (size < 100){printf "hardMasked/%s.hard.fa\n", $1}}' \
| xargs rm
# now get them over to a kluster location
mkdir /san/sanvol1/scratch/bosTau4/hardChunks
cd /san/sanvol1/scratch/bosTau4/hardChunks
# creating 4,000,000 sized chunks, the chroms stay together as
# single pieces. The contigs get grouped together into 4,000,000
# sized fasta files. You don't want to break these things up
# because genscan will be doing its own internal 2.4 million
# window on these pieces, and the gene names are going to be
# constructed from the sequence name in these fasta files.
echo /cluster/data/bosTau4/hardMasked/*.hard.fa | xargs cat \
| faSplit about stdin 4000000 c_
ssh hgwdev
mkdir /cluster/data/bosTau4/bed/genscan
cd /cluster/data/bosTau4/bed/genscan
# Check out hg3rdParty/genscanlinux to get latest genscan:
cvs co hg3rdParty/genscanlinux
# Run on small cluster (more mem than big cluster).
ssh memk
cd /cluster/data/bosTau4/bed/genscan
# Make 3 subdirectories for genscan to put their output files in
mkdir gtf pep subopt
# Generate a list file, genome.list, of all the hard-masked contigs that
# *do not* consist of all-N's (which would cause genscan to blow up)
# Since we split on gaps, we have no chunks like that. You can
# verify with faCount on the chunks.
ls -1Sr /san/sanvol1/scratch/bosTau4/hardChunks/c_*.fa > genome.list
# Create run-time script to operate gsBig in a cluster safe manner
cat << '_EOF_' > runGsBig
#!/bin/csh -fe
set runDir = `pwd`
set srcDir = $1
set inFile = $2
set fileRoot = $inFile:r
mkdir /scratch/tmp/$fileRoot
cp -p $srcDir/$inFile /scratch/tmp/$fileRoot
pushd /scratch/tmp/$fileRoot
/cluster/bin/x86_64/gsBig $inFile $fileRoot.gtf -trans=$fileRoot.pep -subopt=$fileRoot.bed -exe=$runDir/hg3rdParty/genscanlinux/genscan -par=$runDir/hg3rdParty/genscanlinux/HumanIso.smat -tmp=/scratch/tmp -window=2400000
popd
cp -p /scratch/tmp/$fileRoot/$fileRoot.gtf gtf
cp -p /scratch/tmp/$fileRoot/$fileRoot.pep pep
cp -p /scratch/tmp/$fileRoot/$fileRoot.bed subopt
rm -fr /scratch/tmp/$fileRoot
'_EOF_'
# << happy emacs
chmod +x runGsBig
cat << '_EOF_' > template
#LOOP
runGsBig /san/sanvol1/scratch/bosTau4/hardChunks $(file1) {check out line gtf/$(root1).gtf} {check out line pep/$(root1).pep} {check out line subopt/$(root1).bed}
#ENDLOOP
'_EOF_'
# << happy emacs
gensub2 genome.list single template jobList
para create jobList
para try, check, push, check, ...
XXX - running 2008-03-11 11:20
# Completed: 97 of 97 jobs
# CPU time in finished jobs: 93952s 1565.87m 26.10h 1.09d 0.003 y
# IO & Wait Time: 372s 6.19m 0.10h 0.00d 0.000 y
# Average job time: 972s 16.21m 0.27h 0.01d
# Longest running job: 0s 0.00m 0.00h 0.00d
# Longest finished job: 4870s 81.17m 1.35h 0.06d
# Submission to last job: 7902s 131.70m 2.19h 0.09d
# cat and lift the results into single files
ssh kkstore06
cd /cluster/data/bosTau4/bed/genscan
sort -k1,1 -k4.4n gtf/c_*.gtf > genscan.gtf
sort -k1,1 -k2,2n subopt/c_*.bed > genscanSubopt.bed
cat pep/c_*.pep > genscan.pep
# Load into the database as so:
ssh hgwdev
cd /cluster/data/bosTau4/bed/genscan
ldHgGene bosTau4 -gtf genscan genscan.gtf
# Read 49598 transcripts in 346887 lines in 1 files
# 49598 groups 5059 seqs 1 sources 1 feature types
# 49598 gene predictions
hgPepPred bosTau4 generic genscanPep genscan.pep
hgLoadBed bosTau4 genscanSubopt genscanSubopt.bed
# Loaded 528092 elements of size 6
# let's check the numbers
time nice -n +19 featureBits bosTau4 genscan
# 58224594 bases of 2731830700 (2.131%) in intersection
time nice -n +19 featureBits bosTau3 genscan
# 59251085 bases of 2731807384 (2.169%) in intersection
#########################################################################
# CPGISLANDS (DONE - 2008-03-11 - Hiram)
ssh hgwdev
mkdir /cluster/data/bosTau4/bed/cpgIsland
cd /cluster/data/bosTau4/bed/cpgIsland
# Build software from Asif Chinwalla (achinwal@watson.wustl.edu)
cvs co hg3rdParty/cpgIslands
cd hg3rdParty/cpgIslands
make
# There was a problem in here, in both cpg.c and cpg_lh.c:
# warning: conflicting types for built-in function 'malloc'
# comment out the lines where malloc is declared to get this to build
# gcc readseq.c cpg_lh.c -o cpglh.exe
cd ../..
ln -s hg3rdParty/cpgIslands/cpglh.exe .
# There may be warnings about "bad character" for IUPAC ambiguous
# characters like R, S, etc. Ignore the warnings.
mkdir results
echo ../../hardMasked/*.hard.fa | sed -e "s/ /\n/g" | while read F
do
FA=${F/*\/}
C=${FA/.hard.fa/}
echo "./cpglh.exe ${FA} > results/${C}.cpg"
nice -n +19 ./cpglh.exe ${F} > results/${C}.cpg
done > cpglh.out 2>&1 &
# about 5 minutes
# Transform cpglh output to bed +
cat << '_EOF_' > filter.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);
}
'_EOF_'
# << happy emacs
catDir results \
| awk -f filter.awk | sort -k1,1 -k2,2n > cpgIsland.bed
ssh hgwdev
cd /cluster/data/bosTau4/bed/cpgIsland
hgLoadBed bosTau4 cpgIslandExt -tab \
-sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
# Reading cpgIsland.bed
# Loaded 37595 elements of size 10
featureBits bosTau4 cpgIslandExt
# 24202824 bases of 2731830700 (0.886%) in intersection
featureBits bosTau3 cpgIslandExt
# 24374280 bases of 2731807384 (0.892%) in intersection
#########################################################################
# BLASTZ/CHAIN/NET BOSTAU4 (DONE - 2008-03-11 - Hiram)
ssh kkstore02
screen # use a screen to manage this multi-day job
mkdir /cluster/data/hg18/bed/blastzBosTau4.2008-03-11
cd /cluster/data/hg18/bed/
ln -s blastzBosTau4.2008-03-11 blastz.bosTau4
cd blastzBosTau4.2008-03-11
cat << '_EOF_' > DEF
BLASTZ_M=50
# TARGET: Human Hg18
SEQ1_DIR=/cluster/bluearc/scratch/data/hg18/nib
SEQ1_LEN=/cluster/data/hg18/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000
# QUERY: Cow bosTau4
SEQ2_DIR=/san/sanvol1/scratch/bosTau4/bosTau4.2bit
SEQ2_LEN=/cluster/data/bosTau4/chrom.sizes
# Maximum number of scaffolds that can be lumped together
SEQ2_LIMIT=300
SEQ2_CHUNK=20000000
SEQ2_LAP=0
BASE=/cluster/data/hg18/bed/blastzBosTau4.2008-03-11
TMPDIR=/scratch/tmp
'_EOF_'
# << this line keeps emacs coloring happy
time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \
-bigClusterHub=pk -chainMinScore=3000 -chainLinearGap=medium \
-syntenicNet > do.log 2>&1
# real 611m17.901s
cat fb.hg18.chainBosTau4Link.txt
# 1345348636 bases of 2881515245 (46.689%) in intersection
mkdir /cluster/data/bosTau4/bed/blastz.hg18.swap
cd /cluster/data/bosTau4/bed/blastz.hg18.swap
time nice -n +19 doBlastzChainNet.pl \
/cluster/data/hg18/bed/blastzBosTau4.2008-03-11/DEF \
-bigClusterHub=pk -chainMinScore=3000 -chainLinearGap=medium \
-swap -syntenicNet > do.log 2>&1
# broken down during netSynteny.csh due to too many open files on
# a chainSplit
# To get them split, set up a job on memk to do the split
# via chainFilter: chainFilter -t=${T} bigFile.chain > ${T}.chain
# Where the targets T are from the list:
# grep "^chain " bigFile.chain | awk '{print $3}' | sort -u
# see full example below in Platypus Blastz
###########################################################################
# Blastz Platypus ornAna1 (DONE - 2008-03-11,14 - Hiram)
ssh kkstore06
screen # use screen to control this job
mkdir /cluster/data/bosTau4/bed/blastzOrnAna1.2008-03-11
cd /cluster/data/bosTau4/bed/blastzOrnAna1.2008-03-11
cat << '_EOF_' > DEF
# Cow vs. platypus
BLASTZ_Y=3400
BLASTZ_L=6000
BLASTZ_K=2200
BLASTZ_Q=/cluster/data/blastz/HoxD55.q
BLASTZ_M=50
# TARGET: Cow bosTau4
SEQ1_DIR=/san/sanvol1/scratch/bosTau4/bosTau4.2bit
SEQ1_LEN=/cluster/data/bosTau4/chrom.sizes
# Maximum number of scaffolds that can be lumped together
SEQ1_LIMIT=300
SEQ1_CHUNK=20000000
SEQ1_LAP=0
# QUERY: Platypus ornAna1
SEQ2_DIR=/scratch/data/ornAna1/ornAna1.2bit
SEQ2_LEN=/cluster/data/ornAna1/chrom.sizes
SEQ2_CHUNK=20000000
SEQ2_LIMIT=300
SEQ2_LAP=0
BASE=/cluster/data/bosTau4/bed/blastzOrnAna1.2008-03-11
TMPDIR=/scratch/tmp
'_EOF_'
# << happy emacs
time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \
-chainMinScore=5000 -verbose=2 \
-chainLinearGap=loose -bigClusterHub=kk > do.log 2>&1 &
# real 1233m16.397s
# the pk kluster became free as this job had a long time to go.
# Unsure how long it had to go because the batch became confused
# and couldn't calculate eta, it had waiting=0 jobs ? So, chill
# out the batch, allow the running jobs to complete,
# go to pk and get the jobs completed there.
# There was a difficulty with the para.results file. It lost track
# of 1158 jobs, but they were completed and results exist.
# then continuing:
time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \
-chainMinScore=5000 -verbose=2 -smallClusterHub=memk \
-continue=cat -chainLinearGap=loose -bigClusterHub=pk > cat.log 2>&1 &
# real 66m36.465s
cat fb.bosTau4.chainOrnAna1Link.txt
# 201799073 bases of 2731830700 (7.387%) in intersection
mkdir /cluster/data/ornAna1/bed/blastz.bosTau4.swap
cd /cluster/data/ornAna1/bed/blastz.bosTau4.swap
time nice -n +19 doBlastzChainNet.pl \
/cluster/data/bosTau4/bed/blastzOrnAna1.2008-03-11/DEF \
-chainMinScore=5000 -verbose=2 -smallClusterHub=memk \
-swap -chainLinearGap=loose -bigClusterHub=pk > do.log 2>&1 &
# real 106m17.385s
cat fb.ornAna1.chainBosTau4Link.txt
# 187868681 bases of 1842236818 (10.198%) in intersection
# need to split up the chain file to run the reciprocal net
# there are too many chroms to work with chain split
# do it with chainFilter as a kluster job on memk
ssh memk
mkdir /scratch/tmp/ornAna1
cd /scratch/tmp/ornAna1
for R in 0 1 2 3 4 5 6 7
do
echo $R
ssh mkr0u${R} mkdir /scratch/tmp/ornAna1
ssh mkr0u${R} cp -p /cluster/data/bosTau4/bed/blastz.ornAna1/axtChain/bosTau4.ornAna1.all.chain.gz /scratch/tmp/ornAna1
ssh mkr0u${R} gunzip /scratch/tmp/ornAna1/bosTau4.ornAna1.all.chain.gz
done
zcat /cluster/data/bosTau4/bed/blastz.ornAna1/axtChain/bosTau4.ornAna1.all.chain.gz | grep "^chain " | awk '{print $3}' | sort -u > chain.list
cat << '_EOF_' > splitOne.csh
#!/bin/csh -fe
set tmpDir = /scratch/tmp/ornAna1
set T=$1
chainFilter -t=${T} $tmpDir/bosTau4.ornAna1.all.chain > $tmpDir/${T}.chain
'_EOF_'
# << happy emacs
cat << '_EOF_' > template
#LOOP
./splitOne.csh $(file1)
'_EOF_'
# << happy emacs
gensub2 chain.list single template jobList
para create jobList
para try ... check ... push ... time ...
# fetch the results from tmp directories
for R in 0 1 2 3 4 5 6 7
do
ssh mkr0u${R} "(cd /scratch/tmp/ornAna1; cp -p -u chr*.chain /cluster/data/bosTau4/bed/blastz.ornAna1/axtChain/chain)"
done
############################################################################
# BLATSERVERS ENTRY (DONE - 2008-03-12 - 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 ("bosTau4", "blat2", "17778", "1", "0"); \
INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
VALUES ("bosTau4", "blat2", "17779", "0", "1");' \
hgcentraltest
# test it with some sequence
#########################################################################
## 5-Way Multiz (DONE - 2008-03-18 - Hiram)
##
# the all.chain.gz files were split up via kluster jobs on memk
# in order to get mafSynNet files. Example above in ornAna1 blastz
ssh hgwdev
mkdir /cluster/data/bosTau4/bed/multiz5way
cd /cluster/data/bosTau4/bed/multiz5way
# take the 30-way tree from mm9 and eliminate genomes not in
# this alignment
# rearrange to get bosTau4 on the top of the graph
# paste this tree into the on-line phyloGif tool:
# http://genome.ucsc.edu/cgi-bin/phyloGif
# to create the image for the tree diagram
# select the 9 organisms from the 30-way recently done on mouse mm9
/cluster/bin/phast/tree_doctor \
--prune-all-but Human_hg18,Mouse_mm9,Dog_canFam2,Platypus_ornAna1,Cow_bosTau3 \
/cluster/data/mm9/bed/multiz30way/mm9OnTop.fullNames.nh \
| sed -e "s/bosTau3/bosTau4/g" > 5-way.fullNames.nh
# looks something like this:
(((Mouse_mm9:0.325818,Human_hg18:0.126901):0.019763,
(Dog_canFam2:0.156637,Cow_bosTau4:0.163945):0.031326):0.332197,
Platypus_ornAna1:0.488110);
# rearrange to get Cow at the top:
# this leaves us with:
cat << '_EOF_' > bosTau4.5-way.nh
(((Cow_bosTau4:0.163945,Dog_canFam2:0.156637):0.031326,
(Human_hg18:0.126901,Mouse_mm9:0.325818):0.019763):0.332197,
Platypus_ornAna1:0.488110);
'_EOF_'
# << happy emacs
# create a species list from that file:
sed -e 's/[()]//g; s/ /\n/g; s/,/\n/g' bosTau4.5-way.nh \
| sed -e "s/[ \t]*//g; /^[ \t]$/d; /^$/d" | sort -u \
| sed -e "s/.*_//; s/:.*//" | sort > species.list
# verify that has 5 db names in it
# create a stripped down nh file for use in autoMZ run
echo \
`sed 's/[a-zA-Z0-9]*_//g; s/:0.[0-9]*//g; s/[,;]/ /g' bosTau4.5-way.nh \
| sed -e "s/ / /g"` > tree.5.nh
# that looks like, as a single line:
# (((bosTau4 canFam2) (hg18 mm9)) ornAna1)
# verify all blastz's exists
cat << '_EOF_' > listMafs.csh
#!/bin/csh -fe
cd /cluster/data/bosTau4/bed/multiz5way
foreach db (`grep -v bosTau4 species.list`)
set bdir = /cluster/data/bosTau4/bed/blastz.$db
if (-e $bdir/mafRBestNet/chr1.maf.gz) then
echo "$db mafRBestNet"
else if (-e $bdir/mafSynNet/chr1.maf.gz) then
echo "$db mafSynNet"
else if (-e $bdir/mafNet/chr1.maf.gz) then
echo "$db mafNet"
else
echo "$db mafs not found"
endif
end
'_EOF_'
# << happy emacs
chmod +x ./listMafs.csh
# see what it says:
./listMafs.csh
# canFam2 mafSynNet
# hg18 mafSynNet
# mm9 mafSynNet
# ornAna1 mafSynNet
/cluster/bin/phast/all_dists bosTau4.5-way.nh > 5way.distances.txt
grep -i bostau 5way.distances.txt | sort -k3,3n
# Cow_bosTau4 Dog_canFam2 0.320582
# Cow_bosTau4 Human_hg18 0.341935
# Cow_bosTau4 Mouse_mm9 0.540852
# Cow_bosTau4 Platypus_ornAna1 1.015578
# use the calculated
# distances in the table below to order the organisms and check
# the button order on the browser.
# And if you can fill in the table below entirely, you have
# succeeded in finishing all the alignments required.
#
# featureBits chainLink measures
# chainCalJac1Link chain linearGap
# distance on CalJac1 on other minScore
# 1 0.320582 Dog_canFam2 (% 82.846) (% 78.351) 3000 medium
# 2 0.341935 Human_hg18 (% 74.810) (% 77.648) 3000 medium
# 3 0.540852 Mouse_mm9 (% 76.925) (% 74.694) 3000 medium
# 4 1.015578 Platypus_ornAna1 (% 77.296) (% 76.308) 3000 medium
# create a coherent set of all the mafs involved in this run
mkdir mafLinks
cd mafLinks
ln -s ../../blastz.canFam2/mafSynNet ./canFam2
ln -s ../../blastz.hg18/mafSynNet ./hg18
ln -s ../../blastz.mm9/mafSynNet ./mm9
ln -s ../../blastz.ornAna1/mafSynNet ./ornAna1
# check data size:
du -hscL *
# 982M canFam2
# 940M hg18
# 483M mm9
# 79M ornAna1
# 2.5G total
# copy net mafs to cluster-friendly storage
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way/mafLinks
mkdir -p /san/sanvol1/scratch/bosTau4/multiz5way/mafs
rsync -a --progress --copy-links ./ \
/san/sanvol1/scratch/bosTau4/multiz5way/mafs/
# create a run-time list of contigs to operate on, not all contigs
# exist in all alignments, but we want all contig names used in any
# alignment:
cd /san/sanvol1/scratch/bosTau4/multiz5way/mafs
for D in *
do
cd "${D}"
find . -type f
cd ..
done | sed -e "s#^./##; s#.maf##" | sort -u > /tmp/5-way.contig.list
wc -l /tmp/5-way.contig.list
# 2321 /tmp/5-way.contig.list
# ready for the multiz run
ssh pk
mkdir /cluster/data/bosTau4/bed/multiz5way/splitRun
cd /cluster/data/bosTau4/bed/multiz5way/splitRun
scp -p kkstore06:/tmp/5-way.contig.list .
mkdir -p maf run
cd run
mkdir penn
# use latest penn utilities
P=/cluster/bin/penn/multiz.v11.2007-03-19/multiz-tba
cp -p $P/{autoMZ,multiz,maf_project} penn
# set the db and pairs directories here
cat << '_EOF_' > autoMultiz.csh
#!/bin/csh -ef
set db = bosTau4
set c = $1
set result = $2
set resultDir = $result:h
set run = `pwd`
set tmp = /scratch/tmp/$db/multiz.$c
set pairs = /san/sanvol1/scratch/$db/multiz5way/mafs
rm -fr $tmp
mkdir -p $tmp
mkdir -p $resultDir
cp ../../tree.5.nh ../../species.list $tmp
pushd $tmp
foreach s (`grep -v $db species.list`)
set in = $pairs/$s/$c.maf
set out = $db.$s.sing.maf
if (-e $in.gz) then
zcat $in.gz > $out
else if (-e $in) then
cp $in $out
else
echo "##maf version=1 scoring=autoMZ" > $out
endif
end
set path = ($run/penn $path); rehash
$run/penn/autoMZ + T=$tmp E=$db "`cat tree.5.nh`" $db.*.sing.maf $c.maf
popd
cp $tmp/$c.maf $result
rm -fr $tmp
rmdir --ignore-fail-on-non-empty /scratch/tmp/$db
'_EOF_'
# << happy emacs
chmod +x autoMultiz.csh
cat << '_EOF_' > template
#LOOP
./autoMultiz.csh $(root1) {check out line+ /cluster/data/bosTau4/bed/multiz5way/splitRun/maf/$(root1).maf}
#ENDLOOP
'_EOF_'
# << emacs
gensub2 ../5-way.contig.list single template jobList
para create jobList
para try ... check ... push ... etc
# Completed: 2321 of 2321 jobs
# CPU time in finished jobs: 35680s 594.67m 9.91h 0.41d 0.001 y
# IO & Wait Time: 9398s 156.63m 2.61h 0.11d 0.000 y
# Average job time: 19s 0.32m 0.01h 0.00d
# Longest finished job: 1982s 33.03m 0.55h 0.02d
# Submission to last job: 3983s 66.38m 1.11h 0.05d
# put the split maf results back together into a single maf file
# eliminate duplicate comments
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way
grep "^##maf version" splitRun/maf/chr1.maf \
| sort -u > bosTau4.5way.maf
for F in `find ./splitRun/maf -type f -depth`
do
grep -h "^#" "${F}" | egrep -v "maf version=1|eof maf" \
| sed -e "s#/_MZ_[^ ]* # #g; s#__[0-9]##g"
done | sort -u >> bosTau4.5way.maf
for F in `find ./splitRun/maf -type f -depth`
do
grep -v -h "^#" "${F}"
done >> bosTau4.5way.maf
grep "^##eof maf" splitRun/maf/chr1.maf \
| sort -u >> bosTau4.5way.maf
# load tables for a look
ssh hgwdev
mkdir -p /gbdb/bosTau4/multiz5way/maf
ln -s /cluster/data/bosTau4/bed/multiz5way/*.maf \
/gbdb/bosTau4/multiz5way/maf/multiz5way.maf
# this generates an immense multiz5way.tab file in the directory
# where it is running. Best to run this over in scratch.
cd /scratch/tmp
time nice -n +19 hgLoadMaf \
-pathPrefix=/gbdb/bosTau4/multiz5way/maf bosTau4 multiz5way
# real 3m31.593s
# Loaded 5392296 mafs in 1 files from /gbdb/bosTau4/multiz5way/maf
# load summary table
time nice -n +19 cat /gbdb/bosTau4/multiz5way/maf/*.maf \
| hgLoadMafSummary bosTau4 -minSize=30000 -mergeGap=1500 \
-maxSize=200000 multiz5waySummary stdin
# real 3m41.663s
# Created 779088 summary blocks from 10697864 components and 5164426 mafs from stdin
# Gap Annotation
# prepare bed files with gap info
ssh kkstore06
mkdir /cluster/data/bosTau4/bed/multiz5way/anno
cd /cluster/data/bosTau4/bed/multiz5way/anno
mkdir maf run
# these actually already all exist from previous multiple alignments
# remove the echo in front of the twoBitInfo to actually make it work
for DB in `cat ../species.list`
do
CDIR="/cluster/data/${DB}"
if [ ! -f ${CDIR}/${DB}.N.bed ]; then
echo "creating ${DB}.N.bed"
echo twoBitInfo -nBed ${CDIR}/${DB}.2bit ${CDIR}/${DB}.N.bed
else
ls -og ${CDIR}/${DB}.N.bed
fi
done
cd run
rm -f nBeds sizes
for DB in `grep -v bosTau4 ../../species.list`
do
echo "${DB} "
ln -s /cluster/data/${DB}/${DB}.N.bed ${DB}.bed
echo ${DB}.bed >> nBeds
ln -s /cluster/data/${DB}/chrom.sizes ${DB}.len
echo ${DB}.len >> sizes
done
ssh memk
# temporarily copy the bosTau4.5way.maf file onto the memk
# nodes /scratch/data/bosTau4/maf/ directory
for R in 0 1 2 3 4 5 6 7
do
ssh mkr0u${R} rsync -a --progress \
/cluster/data/bosTau4/bed/multiz5way/bosTau4.5way.maf \
/scratch/data/bosTau4/maf/
done
mkdir /cluster/data/bosTau4/bed/multiz5way/anno/splitMaf
# need to split up the single maf file into individual
# per-scaffold maf files to run annotation on
cd /cluster/data/bosTau4/bed/multiz5way/anno/splitMaf
# create bed files to list approximately 1553 scaffolds in
# a single list, approximately 33 lists
cat << '_EOF_' > mkBedLists.pl
#!/usr/bin/env perl
use strict;
use warnings;
my $bedCount = 0;
my $i = 0;
my $bedFile = sprintf("file_%d.bed", $bedCount);
open (BF,">$bedFile") or die "can not write to $bedFile $!";
open (FH,"</cluster/data/bosTau4/chrom.sizes") or
die "can not read /cluster/data/bosTau4/chrom.sizes $!";
while (my $line = <FH>) {
chomp $line;
if ( (($i + 1) % 1553) == 0 ) {
printf "%s\n", $line;
close (BF);
++$bedCount;
$bedFile = sprintf("file_%d.bed", $bedCount);
open (BF,">$bedFile") or die "can not write to $bedFile $!";
}
++$i;
my ($chr, $size) = split('\s+',$line);
printf BF "%s\t0\t%d\t%s\n", $chr, $size, $chr;
}
close (FH);
close (BH);
'_EOF_'
# << happy emacs
chmod +x mkBedLists.pl
./mkBedLists.pl
# now, run a mafsInRegion on each one of those lists
cat << '_EOF_' > runOne
#!/bin/csh -fe
set runDir = "/cluster/data/bosTau4/bed/multiz5way/anno/splitMaf"
set resultDir = $1
set bedFile = $resultDir.bed
mkdir -p $resultDir
mkdir -p /scratch/tmp/bosTau4/$resultDir
pushd /scratch/tmp/bosTau4/$resultDir
mafsInRegion $runDir/$bedFile -outDir . \
/scratch/data/bosTau4/maf/bosTau4.5way.maf
popd
rsync -q -a /scratch/tmp/bosTau4/$resultDir/ ./$resultDir/
rm -fr /scratch/tmp/bosTau4/$resultDir
rmdir --ignore-fail-on-non-empty /scratch/tmp/bosTau4
'_EOF_'
# << happy emacs
chmod +x runOne
cat << '_EOF_' > template
#LOOP
./runOne $(root1)
#ENDLOOP
'_EOF_'
# << happy emacs
ls file*.bed > runList
gensub2 runList single template jobList
para create jobList
para try ... check ... push ... etc
# Completed: 8 of 8 jobs
# CPU time in finished jobs: 1137s 18.96m 0.32h 0.01d 0.000 y
# IO & Wait Time: 552s 9.19m 0.15h 0.01d 0.000 y
# Average job time: 211s 3.52m 0.06h 0.00d
# Longest finished job: 599s 9.98m 0.17h 0.01d
# Submission to last job: 599s 9.98m 0.17h 0.01d
cd /cluster/data/bosTau4/bed/multiz5way/anno/run
cat << '_EOF_' > doAnno.csh
#!/bin/csh -ef
set outDir = ../maf/$2
set result = $3
set input = $1
mkdir -p $outDir
cat $input | \
nice mafAddIRows -nBeds=nBeds stdin /scratch/data/bosTau4/bosTau4.2bit $result
'_EOF_'
# << happy emacs
chmod +x doAnno.csh
cat << '_EOF_' > template
#LOOP
./doAnno.csh $(path1) $(lastDir1) {check out line+ ../maf/$(lastDir1)/$(root1).maf}
#ENDLOOP
'_EOF_'
# << happy emacs
find ../splitMaf -type f -name "*.maf" > maf.list
gensub2 maf.list single template jobList
para create jobList
para try ... check ... push ... etc.
# Completed: 2320 of 2320 jobs
# CPU time in finished jobs: 8405s 140.09m 2.33h 0.10d 0.000 y
# IO & Wait Time: 6934s 115.56m 1.93h 0.08d 0.000 y
# Average job time: 7s 0.11m 0.00h 0.00d
# Longest finished job: 369s 6.15m 0.10h 0.00d
# Submission to last job: 834s 13.90m 0.23h 0.01d
# put the results back together into a single file
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way/anno
grep "^##maf version" maf/file_0/chr1.maf \
| sort -u > bosTau4.anno.5way.maf
find ./maf -type f -depth -name "*.maf" | while read F
do
grep -v -h "^#" "${F}"
done >> bosTau4.anno.5way.maf
echo "##eof maf" >> bosTau4.anno.5way.maf
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/anno
mkdir -p /gbdb/bosTau4/multiz5way/anno
ln -s `pwd`/bosTau4.anno.5way.maf \
/gbdb/bosTau4/multiz5way/anno/multiz5way.maf
# by loading this into the table multiz5way, it will replace the
# previously loaded table with the unannotated mafs
# huge temp files are made, do them on local disk
cd /scratch/tmp
time nice -n +19 hgLoadMaf -pathPrefix=/gbdb/bosTau4/multiz5way/anno \
bosTau4 multiz5way
# Loaded 6711605 mafs in 1 files from /gbdb/bosTau4/multiz5way/anno
# real 4m54.025s
# normally filter this for chrom size > 1,000,000 and only load
# those chroms. But this is a scaffold assembly, load everything:
time nice -n +19 hgLoadMafSummary bosTau4 -minSize=30000 -mergeGap=1500 \
-maxSize=200000 multiz5waySummary \
/gbdb/bosTau4/multiz5way/anno/multiz5way.maf
# Created 779088 summary blocks from 10697864 components
# and 6412983 mafs from /gbdb/bosTau4/multiz5way/anno/multiz5way.maf
# real 4m6.782s
# by loading this into the table multiz5waySummary, it will replace
# the previously loaded table with the unannotated mafs
# remove the multiz5way*.tab files in this /scratch/tmp directory
rm multiz5way*.tab
# And, you can remove the previously loaded non-annotated maf file link:
rm /gbdb/bosTau4/multiz5way/maf/bosTau4.5way.maf
rmdir /gbdb/bosTau4/multiz5way/maf
###########################################################################
## Annotate 5-way multiple alignment with gene annotations
## (DONE - 2008-03-18 - Hiram)
# Gene frames
## given previous survey done for 9-way alignment on Marmoset
## and the appearance of new ensGene tables on everything
# use knownGene for hg18, mm9
# use ensGene for canFam2, ornAna1
# and refGene for bosTau4
ssh hgwdev
mkdir /cluster/data/bosTau4/bed/multiz5way/frames
cd /cluster/data/bosTau4/bed/multiz5way/frames
mkdir genes
# knownGene
for DB in hg18 mm9
do
hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from knownGene" ${DB} \
| genePredSingleCover stdin stdout | gzip -2c \
> /scratch/tmp/${DB}.tmp.gz
mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz
echo "${DB} done"
done
# ensGene
for DB in canFam2 ornAna1
do
hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ensGene" ${DB} \
| genePredSingleCover stdin stdout | gzip -2c \
> /scratch/tmp/${DB}.tmp.gz
mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz
echo "${DB} done"
done
# refGene
for DB in bosTau4
do
hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from refGene" ${DB} \
| genePredSingleCover stdin stdout | gzip -2c \
> /scratch/tmp/${DB}.tmp.gz
mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz
echo "${DB} done"
done
ls -og genes
# -rw-rw-r-- 1 861210 Mar 18 15:40 bosTau4.gp.gz
# -rw-rw-r-- 1 1865308 Mar 18 15:40 canFam2.gp.gz
# -rw-rw-r-- 1 2008806 Mar 18 15:39 hg18.gp.gz
# -rw-rw-r-- 1 1965274 Mar 18 15:39 mm9.gp.gz
# -rw-rw-r-- 1 1347532 Mar 18 15:40 ornAna1.gp.gz
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way/frames
# anything to annotate is in a pair, e.g.: bosTau4 genes/bosTau4.gp.gz
time (cat ../anno/bosTau4.anno.5way.maf | nice -n +19 genePredToMafFrames bosTau4 stdin stdout bosTau4 genes/bosTau4.gp.gz hg18 genes/hg18.gp.gz mm9 genes/mm9.gp.gz canFam2 genes/canFam2.gp.gz ornAna1 genes/ornAna1.gp.gz | gzip > multiz5way.mafFrames.gz) > frames.log 2>&1
# see what it looks like in terms of number of annotations per DB:
zcat multiz5way.mafFrames.gz | cut -f4 | sort | uniq -c | sort -n
# 77526 bosTau4
# 110768 ornAna1
# 221524 mm9
# 230010 hg18
# 243396 canFam2
# load the resulting file
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/frames
time nice -n +19 hgLoadMafFrames bosTau4 multiz5wayFrames \
multiz5way.mafFrames.gz
# real 0m21.968s
# enable the trackDb entries:
# frames multiz5wayFrames
# irows on
#############################################################################
# phastCons 5-way (DONE - 2008-03-19 - Hiram)
# split 5way mafs into 10M chunks and generate sufficient statistics
# files for # phastCons
ssh memk
mkdir /cluster/data/bosTau4/bed/multiz5way/msa.split
cd /cluster/data/bosTau4/bed/multiz5way/msa.split
mkdir -p /san/sanvol1/scratch/bosTau4/multiz5way/cons/ss
cat << '_EOF_' > doSplit.csh
#!/bin/csh -ef
set MAFS = /cluster/data/bosTau4/bed/multiz5way/anno/maf
set WINDOWS = /san/sanvol1/scratch/bosTau4/multiz5way/cons/ss
pushd $WINDOWS
set resultDir = $1
set c = $2
rm -fr $resultDir/$c
mkdir -p $resultDir
twoBitToFa -seq=$c /scratch/data/bosTau4/bosTau4.2bit /scratch/tmp/bosTau4.$c.fa
# need to truncate odd-ball scaffold/chrom names that include dots
# as phastCons utils can't handle them
set TMP = /scratch/tmp/$c.clean.maf.$$
perl -wpe 's/^s ([^.]+\.[^. ]+)\.\S+/s $1/' $MAFS/$resultDir/$c.maf > $TMP
/cluster/bin/phast/$MACHTYPE/msa_split $TMP -i MAF \
-M /scratch/tmp/bosTau4.$c.fa \
-o SS -r $resultDir/$c -w 10000000,0 -I 1000 -B 5000
rm -f /scratch/tmp/bosTau4.$c.fa
rm -f $TMP
popd
mkdir -p $resultDir
date > $resultDir/$c.out
'_EOF_'
# << happy emacs
chmod +x doSplit.csh
cat << '_EOF_' > template
#LOOP
doSplit.csh $(dir1) $(root1) {check out line+ $(dir1)/$(root1).out}
#ENDLOOP
'_EOF_'
# << happy emacs
# create list of maf files:
(cd ../anno/maf; find . -type f) | sed -e "s#^./##" > maf.list
gensub2 maf.list single template jobList
para create jobList
para try ... check ... etc
# Completed: 2320 of 2320 jobs
# CPU time in finished jobs: 1710s 28.50m 0.47h 0.02d 0.000 y
# IO & Wait Time: 6951s 115.85m 1.93h 0.08d 0.000 y
# Average job time: 4s 0.06m 0.00h 0.00d
# Longest finished job: 128s 2.13m 0.04h 0.00d
# Submission to last job: 1048s 17.47m 0.29h 0.01d
# take the cons and noncons trees from the mouse 30-way
# Estimates are not easy to make, probably more correctly,
# take the 30-way .mod file, and re-use it here.
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way
cp -p /cluster/data/mm9/bed/multiz30way/mm9.30way.mod .
# Run phastCons
# This job is I/O intensive in its output files, thus it is all
# working over in /scratch/tmp/
ssh memk
mkdir -p /cluster/data/bosTau4/bed/multiz5way/cons/run.cons
cd /cluster/data/bosTau4/bed/multiz5way/cons/run.cons
# there are going to be several different phastCons runs using
# this same script. They trigger off of the current working directory
# $cwd:t which is the "grp" in this script. It is one of:
# all gliers placentals
# Well, that's what it was when used in the Mm9 30-way,
# in this instance, there is only the directory "all"
cat << '_EOF_' > doPhast.csh
#!/bin/csh -fe
set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast/bin
set subDir = $1
set f = $2
set c = $2:r
set len = $3
set cov = $4
set rho = $5
set grp = $cwd:t
set tmp = /scratch/tmp/$f
set cons = /cluster/data/bosTau4/bed/multiz5way/cons
mkdir -p $tmp
set san = /san/sanvol1/scratch/bosTau4/multiz5way/cons
if (-s $cons/$grp/$grp.non-inf) then
cp -p $cons/$grp/$grp.mod $cons/$grp/$grp.non-inf $tmp
cp -p $san/ss/$subDir/$f.ss $cons/$grp/$grp.mod $cons/$grp/$grp.non-inf $tmp
else
cp -p $cons/$grp/$grp.mod $tmp
cp -p $san/ss/$subDir/$f.ss $cons/$grp/$grp.mod $tmp
endif
pushd $tmp > /dev/null
if (-s $grp.non-inf) then
$PHASTBIN/phastCons $f.ss $grp.mod \
--rho $rho --expected-length $len --target-coverage $cov --quiet \
--not-informative `cat $grp.non-inf` \
--seqname $c --idpref $c --most-conserved $f.bed --score > $f.pp
else
$PHASTBIN/phastCons $f.ss $grp.mod \
--rho $rho --expected-length $len --target-coverage $cov --quiet \
--seqname $c --idpref $c --most-conserved $f.bed --score > $f.pp
endif
popd > /dev/null
mkdir -p $san/$grp/pp/$subDir $san/$grp/bed/$subDir
sleep 4
touch $san/$grp/pp/$subDir $san/$grp/bed/$subDir
rm -f $san/$grp/pp/$subDir/$f.pp
rm -f $san/$grp/bed/$subDir/$f.bed
mv $tmp/$f.pp $san/$grp/pp/$subDir
mv $tmp/$f.bed $san/$grp/bed/$subDir
rm -fr $tmp
'_EOF_'
# << happy emacs
chmod a+x doPhast.csh
# Create parasol batch and run it
pushd /san/sanvol1/scratch/bosTau4/multiz5way/cons
find ./ss -type f -name "*.ss" | sed -e "s#^./##; s/.ss$//" \
> /cluster/data/bosTau4/bed/multiz5way/cons/ss.list
popd
# run for all species
cd ..
mkdir -p all run.cons/all
cd all
/cluster/bin/phast.cz/tree_doctor ../../mm9.30way.mod \
--prune-all-but=bosTau3,hg18,mm9,canFam2,ornAna1 \
| sed -e "s/bosTau3/bosTau4/" > all.mod
cd ../run.cons/all
# root1 == chrom name, file1 == ss file name without .ss suffix
# Create template file for "all" run
cat << '_EOF_' > template
#LOOP
../doPhast.csh $(lastDir1) $(file1) 45 .3 .31 {check out line+ /san/sanvol1/scratch/bosTau4/multiz5way/cons/all/pp/$(lastDir1)/$(file1).pp}
#ENDLOOP
'_EOF_'
# << happy emacs
gensub2 ../../ss.list single template jobList
para create jobList
para try ... check ... push ... etc.
# Completed: 2569 of 2569 jobs
# CPU time in finished jobs: 8636s 143.93m 2.40h 0.10d 0.000 y
# IO & Wait Time: 17371s 289.52m 4.83h 0.20d 0.001 y
# Average job time: 10s 0.17m 0.00h 0.00d
# Longest finished job: 44s 0.73m 0.01h 0.00d
# Submission to last job: 1008s 16.80m 0.28h 0.01d
# create Most Conserved track
ssh kolossus
cd /san/sanvol1/scratch/bosTau4/multiz5way/cons/all
find ./bed -type f -name "chr*.bed" | xargs cat \
| sort -k1,1 -k2,2n | \
awk '{printf "%s\t%d\t%d\tlod=%d\t%s\n", $1, $2, $3, $5, $5;}' | \
/cluster/bin/scripts/lodToBedScore /dev/stdin > mostConserved.bed
# ~ 3 minutes
cp -p mostConserved.bed /cluster/data/bosTau4/bed/multiz5way/cons/all
# load into database
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/cons/all
time nice -n +19 hgLoadBed bosTau4 phastConsElements5way mostConserved.bed
# Loaded 1005876 elements of size 5
# Try for 5% overall cov, and 70% CDS cov
# We don't have any gene tracks to compare CDS coverage
# --rho .31 --expected-length 45 --target-coverage .3
featureBits bosTau4 phastConsElements5way
# 132010504 bases of 2731830700 (4.832%) in intersection
# Create merged posterier probability file and wiggle track data files
# currently doesn't matter where this is performed, the san is the same
# network distance from all machines.
# sort by chromName, chromStart so that items are in numerical order
# for wigEncode
cd /san/sanvol1/scratch/bosTau4/multiz5way/cons/all
mkdir -p phastCons5wayScores
for D in `ls -1d pp/file* | sort -t_ -k2n`
do
TOP=`pwd`
F=${D/pp\/}
out=${TOP}/phastCons5wayScores/${F}.data.gz
echo "${D} > ${F}.data.gz"
cd ${D}
find . -name "*.pp" -type f \
| sed -e "s#^./##; s/chrUn.004./chrUn_004_/; s/-/.-./" \
| sort -t '.' -k1,1 -k3.3n \
| sed -e "s/.-./-/; s/chrUn_004_/chrUn.004./" | xargs cat \
| gzip > ${out}
cd "${TOP}"
done
# copy those files to the downloads area:
# /cluster/data/bosTau4/bed/multiz5way/downloads/phastCons5way/phastConsScores
# for hgdownload downloads
# Create merged posterier probability file and wiggle track data files
# currently doesn't matter where this is performed, the san is the same
# network distance from all machines.
cd /san/sanvol1/scratch/bosTau4/multiz5way/cons/all
ls -1 phastCons5wayScores/*.data.gz | sort -t_ -k2n | xargs zcat \
| wigEncode -noOverlap stdin phastCons5way.wig phastCons5way.wib
# Converted stdin, upper limit 1.00, lower limit 0.00
time nice -n +19 cp -p *.wi? /cluster/data/bosTau4/bed/multiz5way/cons/all
# real 0m40.875s
# Load gbdb and database with wiggle.
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/cons/all
ln -s `pwd`/phastCons5way.wib /gbdb/bosTau4/multiz5way/phastCons5way.wib
time nice -n +19 hgLoadWiggle -pathPrefix=/gbdb/bosTau4/multiz5way bosTau4 \
phastCons5way phastCons5way.wig
# real 1m5.667s
# remove garbage
rm wiggle.tab
# Create histogram to get an overview of all the data
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/cons/all
time nice -n +19 hgWiggle -doHistogram \
-hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \
-db=bosTau4 phastCons5way > histogram.data 2>&1
# real 3m37.316s
# create plot of histogram:
cat << '_EOF_' | gnuplot > histo.png
set terminal png small color x000000 xffffff xc000ff x66ff66 xffff00 x00ffff
set size 1.4, 0.8
set key left box
set grid noxtics
set grid ytics
set title " Cow BosTau4 Histogram phastCons5way track"
set xlabel " phastCons5way score"
set ylabel " Relative Frequency"
set y2label " Cumulative Relative Frequency (CRF)"
set y2range [0:1]
set y2tics
set yrange [0:0.02]
plot "histogram.data" using 2:5 title " RelFreq" with impulses, \
"histogram.data" using 2:7 axes x1y2 title " CRF" with lines
'_EOF_'
# << happy emacs
display histo.png &
# These trackDb entries turn on the wiggle phastCons data track:
# type wigMaf 0.0 1.0
# maxHeightPixels 100:40:11
# wiggle phastCons5way
# spanList 1
- # autoScaleDefault Off
+ # autoScale Off
# windowingFunction mean
# pairwiseHeight 12
# yLineOnOff Off
#############################################################################
# Downloads (DONE - 2008-01-11 - Hiram)
# Let's see if the downloads will work
ssh hgwdev
/cluster/data/bosTau4
# expecting to find repeat masker .out file here:
ln -s bed/RepeatMasker/bosTau4.fa.out .
time nice -n +19 /cluster/bin/scripts/makeDownloads.pl \
-workhorse=hgwdev bosTau4 > jkStuff/downloads.log 2>&1
# real 24m3.210s
# failed making upstream sequences:
# featureBits bosTau4 mgcGenes:upstream:1000 -fa=stdout
# setpriority: Permission denied.
# the 'nice' from my bash shell causes trouble inside the csh
# script which uses nice. Finish off the install step manually
# with the mgcGenes upstreams ...
#############################################################################
# PushQ entries (DONE - 2008-01-11 - Hiram)
ssh hgwdev
/cluster/data/bosTau4
/cluster/bin/scripts/makePushQSql.pl bosTau4 > jkStuff/pushQ.sql
# output warnings:
# bosTau4 does not have seq
# bosTau4 does not have gbMiscDiff
# Could not tell (from trackDb, all.joiner and hardcoded lists of supporting
# and genbank tables) which tracks to assign these tables to:
# genscanPep
#############################################################################
# create download files (DONE - 2008-03-19 - Hiram)
ssh hgwdev
cd /cluster/data/bosTau4
ln -s /cluster/data/bosTau4/bed/repeatMasker/bosTau4.fa.out .
makeDownloads.pl bosTau4 > makeDownloads.log 2>&1
# *EDIT* the README files and ensure they are correct
#############################################################################
# PushQ entries (DONE - 2008-03-19 - Hiram)
ssh hgwdev
/cluster/data/bosTau4
/cluster/bin/scripts/makePushQSql.pl bosTau4 > jkStuff/pushQ.sql
# output warnings:
# hgwdev does not have /usr/local/apache/htdocs/goldenPath/bosTau4/liftOver/bosTau4ToBosTau*
# bosTau4 does not have seq
# Could not tell (from trackDb, all.joiner and hardcoded lists of supporting
# and genbank tables) which tracks to assign these tables to:
# genscanPep
# looks like there should be a bosTau3 to bosTau4 liftOver run
###########################################################################
# HUMAN (hg18) PROTEINS TRACK (DONE braney 2008-03-28)
ssh kkstore06
# bash if not using bash shell already
mkdir /cluster/data/bosTau4/blastDb
cd /cluster/data/bosTau4
grep -v chrUn chrom.sizes | awk '{print $1}' > chr.lst
for i in `cat chr.lst`;
do twoBitToFa bosTau4.unmasked.2bit -seq=$i stdout;
done > temp.fa
faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
grep chrUn chrom.sizes | awk '{print $1}' > chr.lst
for i in `cat chr.lst`;
do twoBitToFa bosTau4.unmasked.2bit -seq=$i stdout;
done > temp.fa
faSplit sequence temp.fa 150 blastDb/y
rm temp.fa chr.lst
cd blastDb
for i in *.fa
do
/cluster/bluearc/blast229/formatdb -i $i -p F
done
rm *.fa
ls *.nsq | wc -l
# 3440
mkdir -p /san/sanvol1/scratch/bosTau4/blastDb
cd /cluster/data/bosTau4/blastDb
for i in nhr nin nsq;
do
echo $i
cp *.$i /san/sanvol1/scratch/bosTau4/blastDb
done
mkdir -p /cluster/data/bosTau4/bed/tblastn.hg18KG
cd /cluster/data/bosTau4/bed/tblastn.hg18KG
echo /san/sanvol1/scratch/bosTau4/blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst
wc -l query.lst
# 3440 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/3440) = 360.973943
mkdir -p /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/kgfa
split -l 361 /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/kgfa/kg
ln -s /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/kgfa kgfa
cd kgfa
for i in *; do
nice pslxToFa $i $i.fa;
rm $i;
done
cd ..
ls -1S kgfa/*.fa > kg.lst
mkdir -p /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/blastOut
ln -s /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/blastOut
for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done
tcsh
cd /cluster/data/bosTau4/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=/cluster/bluearc/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 /cluster/bluearc/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/bosTau4/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
gensub2 query.lst kg.lst blastGsub blastSpec
exit
ssh pk
cd /cluster/data/bosTau4/bed/tblastn.hg18KG
para create blastSpec
# para try, check, push, check etc.
para time
Completed: 350880 of 350880 jobs
CPU time in finished jobs: 27082816s 451380.27m 7523.00h 313.46d 0.859 y
IO & Wait Time: 2334990s 38916.50m 648.61h 27.03d 0.074 y
Average job time: 84s 1.40m 0.02h 0.00d
Longest finished job: 578s 9.63m 0.16h 0.01d
Submission to last job: 96125s 1602.08m 26.70h 1.11d
ssh kkstore06
cd /cluster/data/bosTau4/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 /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/blastOut/c.`basename $1`.psl)
'_EOF_'
chmod +x chainOne
ls -1dS /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/blastOut/kg?? > chain.lst
gensub2 chain.lst single chainGsub chainSpec
# do the cluster run for chaining
ssh pk
cd /cluster/data/bosTau4/bed/tblastn.hg18KG/chainRun
para create chainSpec
para maxNode 30
para try, check, push, check etc.
# Completed: 99 of 102 jobs
# Crashed: 3 jobs
# CPU time in finished jobs: 113248s 1887.47m 31.46h 1.31d 0.004 y
# IO & Wait Time: 86043s 1434.04m 23.90h 1.00d 0.003 y
# Average job time: 2013s 33.55m 0.56h 0.02d
# Longest finished job: 6139s 102.32m 1.71h 0.07d
# Submission to last job: 10416s 173.60m 2.89h 0.12d
# ran three crashed jobs on kolossus
ssh kkstore06
cd /cluster/data/bosTau4/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 -T /tmp -k 14,14 -k 16,16n -k 17,17n u.*.psl m60* | uniq > /cluster/data/bosTau4/bed/tblastn.hg18KG/blastHg18KG.psl
cd ..
pslCheck blastHg18KG.psl
# load table
ssh hgwdev
cd /cluster/data/bosTau4/bed/tblastn.hg18KG
hgLoadPsl bosTau4 blastHg18KG.psl
# check coverage
featureBits bosTau4 blastHg18KG
# 40254923 bases of 2731830700 (1.474%) in intersection
featureBits bosTau4 refGene:cds blastHg18KG -enrichment
# refGene:cds 0.429%, blastHg18KG 1.474%, both 0.379%, cover 88.39%, enrich 59.98x
ssh kkstore06
rm -rf /cluster/data/bosTau4/bed/tblastn.hg18KG/blastOut
rm -rf /cluster/bluearc/bosTau4/bed/tblastn.hg18KG/blastOut
#end tblastn
##########################################################################
# Create 5-way downloads (DONE - 2008-03-28 - Hiram)
ssh hgwdev
mkdir -p /cluster/data/bosTau4/bed/multiz5way/downloads/phastCons5way
cd /cluster/data/bosTau4/bed/multiz5way/downloads/phastCons5way
cp -p \
/san/sanvol1/scratch/bosTau4/multiz5way/cons/all/phastCons5wayScores/* .
ln -s ../../cons/all/all.mod ./5way.mod
cp /cluster/data/calJac1/bed/multiz9way/downloads/phastCons9way/README.txt .
# edit that README.txt to be correct for this 5-way alignment
cd ..
mkdir multiz5way
cd multiz5way
cp -p /cluster/data/calJac1/bed/multiz9way/downloads/multiz9way/README.txt .
# edit that README.txt to be correct for this 5-way alignment
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way/downloads/multiz5way
ln -s ../../bosTau4.5-way.nh 5way.nh
time gzip -c ../../anno/bosTau4.anno.5way.maf > bosTau4.5way.maf.gz
# real 34m59.295s
ssh hgwdev
cd /cluster/data/bosTau4/bed/multiz5way/downloads/multiz5way
# creating upstream files from refGene, bash script:
cat << '_EOF_' > mkUpstream.sh
#!/bin/bash
DB=bosTau4
GENE=refGene
NWAY=multiz5way
export DB GENE
for S in 1000 2000 5000
do
echo "making upstream${S}.maf"
featureBits ${DB} ${GENE}:upstream:${S} -fa=/dev/null -bed=stdout \
| perl -wpe 's/_up[^\t]+/\t0/' | sort -k1,1 -k2,2n \
| $HOME/kent/src/hg/ratStuff/mafFrags/mafFrags ${DB} ${NWAY} \
stdin stdout \
-orgs=/cluster/data/${DB}/bed/${NWAY}/species.list \
| gzip -c > upstream${S}.maf.gz
echo "done upstream${S}.maf.gz"
done
'_EOF_'
# << happy emacs
chmod +x ./mkUpstream.sh
time nice -n +19 ./mkUpstream.sh
-rw-rw-r-- 1 9883443 Mar 28 13:02 upstream1000.maf.gz
-rw-rw-r-- 1 17938570 Mar 28 13:06 upstream2000.maf.gz
-rw-rw-r-- 1 40384656 Mar 28 13:10 upstream5000.maf.gz
#
# check the names in these upstream files to ensure sanity:
zcat upstream1000.maf.gz | grep "^s " | awk '{print $2}' \
| sort | uniq -c | sort -rn | less
# should be a list of the other 4 species with a high count,
# then refGene names, e.g.:
# 8806 ornAna1
# 8806 mm9
# 8806 hg18
# 8806 canFam2
# 7 NM_001077006
# 3 NM_001113231
# 3 NM_001105381
# 3 NM_001102527
# 3 NM_001102322
# ...
ssh kkstore06
cd /cluster/data/bosTau4/bed/multiz5way/downloads/multiz5way
md5sum *.maf.gz > md5sum.txt
cd ../phastCons5way
md5sum *.data.gz *.mod > md5sum.txt
ssh hgwdev
mkdir /usr/local/apache/htdocs/goldenPath/bosTau4/multiz5way
mkdir /usr/local/apache/htdocs/goldenPath/bosTau4/phastCons5way
cd /cluster/data/bosTau4/bed/multiz5way/downloads/multiz5way
ln -s `pwd`/* /usr/local/apache/htdocs/goldenPath/bosTau4/multiz5way
cd ../phastCons5way
ln -s `pwd`/* /usr/local/apache/htdocs/goldenPath/bosTau4/phastCons5way
# if your ln -s `pwd`/* made extra links to files you don't want there,
# check the goldenPath locations and remove those extra links
#############################################################################
# N-SCAN gene predictions (nscanGene) - (2008-04-21 markd)
# obtained NSCAN predictions from michael brent's group
# at WUSTL
cd /cluster/data/bosTau4/bed/nscan/
wget http://mblab.wustl.edu/predictions/Cow/bosTau4/bosTau4.gtf
wget http://mblab.wustl.edu/predictions/Cow/bosTau4/bosTau4.prot.fa
wget http://mblab.wustl.edu/predictions/Cow/bosTau4/readme.html
bzip2 bosTau4.*
chmod a-w *
# load track
gtfToGenePred -genePredExt bosTau4.gtf.bz2 stdout | hgLoadGenePred -bin -genePredExt bosTau4 nscanGene stdin
hgPepPred bosTau4 generic nscanPep bosTau4.prot.fa.bz2
rm *.tab
# update trackDb; need a bosTau4-specific page to describe informants
cow/bosTau4/nscanGene.html (copy from readme.html)
cow/bosTau4/trackDb.ra
# set search regex to
termRegex chr[0-9a-zA-Z_]+\.[0-9]+\.[0-9]+(\.[0-9]+)*
#############################################################################
############################################################################
# TRANSMAP vertebrate.2008-05-20 build (2008-05-24 markd)
vertebrate-wide transMap alignments were built Tracks are created and loaded
by a single Makefile. This is available from:
svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2008-05-20
see doc/builds.txt for specific details.
############################################################################
# SELF BLASTZ (DONE - 2008-06-27 - Hiram)
# do not want the noise from the contigs
ssh kkstore06
screen # use a screen to manage this multi-day job
mkdir /cluster/data/bosTau4/noContigs
cd /cluster/data/bosTau4/noContigs
cut -f1 ../chrom.sizes | grep -v chrUn | while read C
do
echo "twoBitToFa -seq=${C} ../bosTau4.2bit ${C}.fa"
twoBitToFa -seq=${C} ../bosTau4.2bit ${C}.fa
done
faToTwoBit chr*.fa bosTau4.noContigs.2bit
twoBitToFa bosTau4.noContigs.2bit stdout | faSize stdin
# 2634429662 bases (167456864 N's 2466972798 real
# 1326108891 upper 1140863907 lower) in 31 sequences in 1 files
# %43.31 masked total, %46.25 masked real
grep -v chrUn ../chrom.sizes | ave -col=2 stdin
# count 31
# total 2634429662.000000
twoBitInfo bosTau4.noContigs.2bit stdout | sort -k2nr \
> bosTau4.noContigs.chrom.sizes
cp -p bosTau4.noContigs.2bit /cluster/bluearc/scratch/data/bosTau4
cp -p bosTau4.noContigs.chrom.sizes /cluster/bluearc/scratch/data/bosTau4
mkdir /cluster/data/bosTau4/bed/blastzSelf.2008-08-27
cd /cluster/data/bosTau4/bed
ln -s blastzSelf.2008-08-27 blastzSelf
cd blastzSelf.2008-08-27
cat << '_EOF_' > DEF
BLASTZ_M=400
# TARGET: Cow bosTau4
SEQ1_DIR=/cluster/bluearc/scratch/data/bosTau4/bosTau4.noContigs.2bit
SEQ1_LEN=/cluster/bluearc/scratch/data/bosTau4/bosTau4.noContigs.chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000
# QUERY: Cow bosTau4
SEQ2_DIR=/cluster/bluearc/scratch/data/bosTau4/bosTau4.noContigs.2bit
SEQ2_LEN=/cluster/bluearc/scratch/data/bosTau4/bosTau4.noContigs.chrom.sizes
SEQ2_CHUNK=10000000
SEQ2_LAP=0
BASE=/cluster/data/bosTau4/bed/blastzSelf.2008-08-27
TMPDIR=/scratch/tmp
'_EOF_'
# << this line keeps emacs coloring happy
cd /cluster/data/bosTau4/bed/blastzSelf.2008-08-27
time /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
-chainMinScore=10000 -chainLinearGap=medium -bigClusterHub=pk \
`pwd`/DEF > blastz.out 2>&1 &
# real 3537m49.719s
# had to fix a slight problem in the make downloads step, then:
time /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
-chainMinScore=10000 -chainLinearGap=medium -bigClusterHub=pk \
-continue=cleanup `pwd`/DEF > cleanup.out 2>&1 &
############################################################################
############################################################################
# TRANSMAP vertebrate.2008-06-07 build (2008-06-30 markd)
vertebrate-wide transMap alignments were built Tracks are created and loaded
by a single Makefile. This is available from:
svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2008-06-30
see doc/builds.txt for specific details.
############################################################################
################################################
# AUTOMATE UPSTREAM FILE CREATION (2008-10-15 markd)
update genbank.conf:
bosTau4.upstreamGeneTbl = mgcGenes
bosTau4.upstreamMaf = multiz5way /hive/data/genomes/bosTau4/bed/multiz5way/species.list
############################################################################
# TRANSMAP vertebrate.2009-07-01 build (2009-07-21 markd)
vertebrate-wide transMap alignments were built Tracks are created and loaded
by a single Makefile. This is available from:
svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2009-07-01
see doc/builds.txt for specific details.
############################################################################
############################################################################
# TRANSMAP vertebrate.2009-09-13 build (2009-09-20 markd)
vertebrate-wide transMap alignments were built Tracks are created and loaded
by a single Makefile. This is available from:
svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2009-09-13
see doc/builds.txt for specific details.
############################################################################