src/hg/makeDb/doc/felCatV17e.txt 1.13
1.13 2010/05/07 17:20:37 chinhli
Complete multiz6way
Index: src/hg/makeDb/doc/felCatV17e.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/felCatV17e.txt,v
retrieving revision 1.12
retrieving revision 1.13
diff -b -B -U 1000000 -r1.12 -r1.13
--- src/hg/makeDb/doc/felCatV17e.txt 5 May 2010 14:56:10 -0000 1.12
+++ src/hg/makeDb/doc/felCatV17e.txt 7 May 2010 17:20:37 -0000 1.13
@@ -1,1180 +1,1300 @@
# for emacs: -*- mode: sh; -*-
# $Id$
# Felis Catus (domestic cat) -- NHGRI/GTB V17e/felCatV17e (2009-05-25)
# file template copied from calJac3.txt
# Felis catus (Project ID: 32759) by NHGRI/Genome Technology Branch [Draft
# assembly] sequence:
# ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/
# Felis_catus/catChrV17e
# Felis catus
##########################################################################
# Download sequence (DONE - 2010-03-03 Chin)
mkdir /hive/data/genomes/felCatV17e
cd /hive/data/genomes/felCatV17e
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_mammals/Felis_catus/catChrV17e/*"
# FINISHED --09:05:15--
# Downloaded: 151 files, 1.3G in 7m 42s (2.98 MB/s)
# Read ASSEMBLY_INFO
mkdir ucscChr
# stay at genbank directory
# fixup the accession names to become UCSC chrom names
S=Primary_Assembly/assembled_chromosomes
cut -f1 ${S}/chr2acc | while read C
do
ACC=`grep "${C}" ${S}/chr2acc | cut -f2`
echo "${ACC} -> chr${C}"
zcat ${S}/AGP/chr${C}.agp.gz \
| sed -e "s/^${ACC}/chr${C}/" | gzip > ucscChr/chr${C}.agp.gz
done
# use component instead
S=Primary_Assembly/assembled_chromosomes
cut -f1 ${S}/chr2acc | while read C
do
ACC=`grep "${C}" ${S}/chr2acc | cut -f2`
echo "${ACC} -> chr${C}"
zcat ${S}/AGP/chr${C}.comp.agp.gz \
| sed -e "s/^${ACC}/chr${C}/" | gzip > ucscChr/chr${C}.comp.agp.gz
done
S=Primary_Assembly/assembled_chromosomes
cut -f1 ${S}/chr2acc | while read C
do
ACC=`grep "${C}" ${S}/chr2acc | cut -f2`
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
# 2872644707 bases (1165972091 N's 1706672616 real 1706672616 upper 0
# lower) in 19 sequences in 19 files
faSize ucscChr/chr*.fa.gz
# 2872644707 bases (1165972091 N's 1706672616 real 1706672616 upper 0
# lower) in 19 sequences in 19 files
# For unplaced scalfolds, named them as chrUn_xxxxxxxx
# where xxxxxx is the original access id as: chrUn_ACBE01000381.1
# and put it into chrUn.* files
# copy all the comment lines (start with #)
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_/" >> ucscChr/chrUn.agp
gzip ucscChr/chrUn.agp &
zcat ${S}/FASTA/unplaced.scaf.fa.gz \
| sed -e "s#^>.*|gb|#>chrUn_#; s#|.*##" | gzip > ucscChr/chrUn.fa.gz
# about 104034 sequences in the unplaced
zcat ucscChr/chrUn.fa.gz | grep "^>" | wc
# 104034 104034 2275648
# Check them with faSize
faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz
# 287642232 bases (3696852 N's 283945380 real 283945380 upper 0
# lower) in 104034 sequences in 1 files
faSize ucscChr/chrUn.fa.gz
# 287642232 bases (3696852 N's 283945380 real 283945380 upper 0
# lower) in 104034 sequences in 1 files
##########################################################################
# Initial genome build (Done - 2010-03-03 -Chin)
# decide the right order by
hgsql -e "select orderKey,name from dbDb order by orderKey;" \
hgcentraltest | less
# 168 cavPor3
# 169 cavPor2
# 189 oryCun2
# 190 oryCun1
# 191 ochPri2
# 200 eriEur1
# 217 felCat3
cd /hive/data/genomes/felCatV17e
cat << '_EOF_' > felCatV17e.config.ra
# Config parameters for makeGenomeDb.pl:
db felCatV17e
clade mammal
genomeCladePriority 16
scientificName Felis catus
commonName Cat
assemblyDate Dec. 2008
assemblyLabel NHGRI/Genome Technology Branch (NCBI project 10703, accession ACBE0100000)
assemblyShortLabel NHGRI/GTB V17e
orderKey 216
mitoAcc NC_001700
fastaFiles /hive/data/genomes/felCatV17e/genbank/ucscChr/chr*.fa.gz
agpFiles /hive/data/genomes/felCatV17e/genbank/ucscChr/chr*.agp.gz
# qualFiles none
dbDbSpeciesDir cat
taxId 9685
'_EOF_'
makeGenomeDb.pl -stop seq felCatV17e.config.ra > seq.log 2>&1 &
time makeGenomeDb.pl -continue=agp -stop=agp felCatV17e.config.ra > agp.log 2>&1 &
# real 0m50.486s
time makeGenomeDb.pl -continue=db -stop=db felCatV17e.config.ra > db.log 2>&1 &
#real 7m50.591s
time makeGenomeDb.pl -continue=dbDb -stop=dbDb felCatV17e.config.ra > dbDb.log 2>&1 &
time makeGenomeDb.pl -continue=trackDb felCatV17e.config.ra > trackDb.log 2>&1 &
#### 03/15 fixing the agp using the component one
# manually run hgLoadGapGl (as in jkStuff/makeDb.csh)
# without rerun the makeGenomeDb -continue=gap step
zcat genbank/ucscChr/chr*.comp.agp.gz | grep -v '^#' \
| sort -k1,1 -k2n,2n > felCatV17e.comp.agp
zcat genbank/ucscChr/chrUn.agp.gz | grep -v '^#' >> felCatV17e.comp.agp
grep chrM felCatV17e.agp >> felCatV17e.comp.agp
checkAgpAndFa felCatV17e.comp.agp felCatV17e.unmasked.2bit
# loop: chrM, dnaOffset=0, seqSize=17009
# chrM 1 17009 230 F NC_001700 1 17009 +
# agpFrag->chromStart: 0, agpFrag->chromEnd: 17009, dnaOffset: 0
# FASTA sequence entry
# valid Fasta file entry
# all AGP and FASTA entries agree - both files are valid
dbSnoop felCatV17e before.txt
hgGoldGapGl -noGl felCatV17e felCatV17e.comp.agp
dbSnoop felCatV17e after.txt
diff before.txt after.txt # to see the change of rows in gold and gap tables
# > 35,861,584 gold 506,133 24,249,424 11,612,160
# > 28,192,744 gap 500,507 22,068,200 6,124,544
# 26,27d27
# < 858,728 gold 9,961 539,240 319,488
# < 345,208 gap 4,335 236,664 108,544
featureBits -countGaps felCatV17e gap
# 1169658726 bases of 3160303948 (37.011%) in intersection
hgsql -e "select type from gap;" felCatV17e | sort | uniq -c
# 381 contig
# 500126 fragment
# 1 type
mv felCatV17e.comp.agp felCatV17e.agp
##########################################################################
# fixed up the html files (working - 2010-04-06 Chin)
# per instructions at end of makeGenomeDb.log, edit the html file
# and checkin the *.ra and *.html files.
cd /cluster/home/chinhli/kent/src/hg/makeDb/trackDb/cat/felCatV17e
##########################################################################
# running repeat masker (DONE - 2010-03-05 - Chin)
mkdir /hive/data/genomes/felCatV17e/bed/repeatMasker
cd /hive/data/genomes/felCatV17e/bed/repeatMasker
time doRepeatMasker.pl -buildDir=`pwd` -noSplit \
-bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \
-smallClusterHub=encodek felCatV17e > do.log 2>&1 &
# real 530m25.613s
cat faSize.rmsk.txt
# 3160303948 bases (1169668943 N's 1990635005 real 1183631401
# upper 807003604 lower) in 104054 sequences in 1 files
# %25.54 masked total, %40.54 masked real
##########################################################################
# running simple repeat (DONE - 2010-03-16 Chin)
mkdir /hive/data/genomes/felCatV17e/bed/simpleRepeat
cd /hive/data/genomes/felCatV17e/bed/simpleRepeat
time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \
-dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=encodek \
felCatV17e > do.log 2>&1 &
# reall 149m19.948s
cat fb.simpleRepeat
# 49107705 bases of 2690950627 (1.825%) in intersection
cd /hive/data/genomes/felCatV17e
twoBitMask felCatV17e.rmsk.2bit \
-add bed/simpleRepeat/trfMask.bed felCatV17e.2bit
# you can safely ignore the warning about fields >= 13
twoBitToFa felCatV17e.2bit stdout | faSize stdin > faSize.felCatV17e.2bit.txt
cat faSize.felCatV17e.2bit.txt
# 3160303948 bases (1169668943 N's 1990635005 real 1182766163 upper
# 807868842 lower) in 104054 sequences in 1 files
# %25.56 masked total, %40.58 masked real
# double check with featureBits
featureBits -countGaps felCatV17e gap
# 1169668943 bases of 3160303948 (37.011%) in intersection
rm /gbdb/felCatV17e/felCatV17e.2bit
ln -s `pwd`/felCatV17e.2bit /gbdb/felCatV17e/felCatV17e.2bit
########################################################################
# Marking *all* gaps - they are not all in the AGP file
# (DONE - 2010-03-16 - Chin)
mkdir /hive/data/genomes/felCatV17e/bed/allGaps
cd /hive/data/genomes/felCatV17e/bed/allGaps
time nice -n +19 findMotif -motif=gattaca -verbose=4 \
-strand=+ ../../felCatV17e.rmsk.2bit > findMotif.txt 2>&1
# real 1m45.351s
grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed
featureBits felCatV17e -not gap -bed=notGap.bed
# 1990645222 bases of 1990645222 (100.000%) in intersection
featureBits felCatV17e allGaps.bed notGap.bed -bed=new.gaps.bed
# 10217 bases of 1990645222 (0.001%) in intersection
# (real time 7+ hours)
# what is the last index in the existing gap table:
hgsql -N -e "select ix from gap;" felCatV17e | sort -n | tail -1
# 101050
cat << '_EOF_' > mkGap.pl
#!/usr/bin/env perl
use strict;
use warnings;
my $ix=`hgsql -N -e "select ix from gap;" felCatV17e | sort -n | tail -1`;
chomp $ix;
open (FH,"<new.gaps.bed") or die "can not read new.gaps.bed";
while (my $line = <FH>) {
my ($chrom, $chromStart, $chromEnd, $rest) = split('\s+', $line);
++$ix;
printf "%s\t%d\t%d\t%d\tN\t%d\tother\tyes\n", $chrom, $chromStart,
$chromEnd, $ix, $chromEnd-$chromStart;
}
close (FH);
'_EOF_'
# << happy emacs
chmod +x ./mkGap.pl
./mkGap.pl > other.gap
hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \
-noLoad felCatV17e otherGap other.gap
# Reading other.gap
# Loaded 1294 elements of size 8
# Sorted
# Saving bed.tab
# No load option selected, see file: bed.tab
wc -l bed.tab
# 1294 bed.tab
# starting with this many
hgsql -e "select count(*) from gap;" felCatV17e
# 500507
hgsql felCatV17e -e 'load data local infile "bed.tab" into table gap;'
# result count:
hgsql -e "select count(*) from gap;" felCatV17e
# 501801
# == 500507 + 1294
#########################################################################
# MAKE 11.OOC FILES FOR BLAT (DONE - 2010-02-03 - Chin)
cat faSize.felCatV17e.2bit.txt
# 3160303948 bases (1169668943 N's 1990635005 real 1182766163 upper
# 807868842 lower) in 104054 sequences in 1 files
# %25.56 masked total, %40.58 masked real
# numerator is felCatV17e gapless bases "real" as reported by faSize
# denominator is hg17 gapless bases as reported by featureBits,
# 1024 is threshold used for human -repMatch:
calc \( 1990635005 / 2897310462 \) \* 1024
# ( 1990635005 / 2897310462 ) * 1024 = 703.552578
# ==> use -repMatch=700 according to size scaled down from 1024 for human.
# and rounded down to nearest 50
cd /hive/data/genomes/felCatV17e
blat felCatV17e.2bit /dev/null /dev/null -tileSize=11 \
-makeOoc=jkStuff/felCatV17e.11.ooc -repMatch=700
# Wrote 23033 overused 11-mers to jkStuff/felCatV17e.11.ooc
mkdir /hive/data/staging/data/felCatV17e
cp -p felCatV17e.2bit chrom.sizes jkStuff/felCatV17e.11.ooc \
/hive/data/staging/data/felCatV17e
gapToLift -bedFile=jkStuff/nonBridgedGaps.bed felCatV17e \
jkStuff/felCatV17e.nonBridged.lft
# Ask the admin to copy the
# /hive/data/staging/data/felCatV17e directory
# to cluster nodes: /scratch/data/felCatV17e
##########################################################################
# BLATSERVERS ENTRY (DONE - 2009-12-23 - Chin)
# After getting a blat server assigned by the Blat Server Gods,
#
hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
VALUES ("felCatV17e", "blat4", "17778", "1", "0"); \
INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
VALUES ("felCatV17e", "blat4", "17779", "0", "1");' \
hgcentraltest
# test it with some sequence
############################################################################
# chrA2:64650765-64656835
# reset position to RHO location as found from blat of hg19 RHO gene
hgsql -e \
'update dbDb set defaultPos="chrA2:64650765-64656835" where
name="felCatV17e";' \
hgcentraltest
############################################################################
# genbank run DONE - 2010-03-16 - Chin )
ssh hgwdev
cd $HOME/kent/src/hg/makeDb/genbank
# edit etc/genbank.conf to add this section just before calJac1:
# Marmoset
felCatV17e.serverGenome = /hive/data/genomes/felCatV17e/felCatV17e.2bit
felCatV17e.clusterGenome = /scratch/data/felCatV17e/felCatV17e.2bit
felCatV17e.ooc = /scratch/data/felCatV17e/felCatV17e.11.ooc
felCatV17e.lift = no
felCatV17e.perChromTables = no
felCatV17e.refseq.mrna.native.pslCDnaFilter = ${ordered.refseq.mrna.native.pslCDnaFilter}
felCatV17e.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter}
felCatV17e.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter}
felCatV17e.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter}
felCatV17e.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter}
felCatV17e.genbank.est.xeno.pslCDnaFilter = ${ordered.genbank.est.xeno.pslCDnaFilter}
felCatV17e.downloadDir = felCatV17e
felCatV17e.refseq.mrna.native.load = yes
felCatV17e.refseq.mrna.xeno.load = yes
felCatV17e.refseq.mrna.xeno.loadDesc = yes
cvs ci -m "adding felCatV17e" 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 felCatV17e &
# tail var/build/logs/2010.03.17-10:43:18.felCatV17e.initalign.log
# gbAlignInstall: complete: real=3.63
# genbank 2010.03.17-16:33:21 felCatV17e.initalign: finish
ssh hgwdev
cd /cluster/data/genbank
time ./bin/gbDbLoadStep -drop -initialLoad felCatV17e &
# logFile:
# tail var/dbload/hgwdev/logs/2010.03.18-15:03:58.dbload.log
# real 31m56.710s
# enable daily alignment and update of hgwdev
cd ~/kent/src/hg/makeDb/genbank
cvsup
# add felCatV17e to:
etc/align.dbs
etc/hgwdev.dbs
cvs ci -m "Adding felCatV17e - Cat - Felis catus" \
etc/align.dbs etc/hgwdev.dbs
make etc-update
# DONE - 2010-03-17 - Chin
############################################################################
# running cpgIsland business (DONE - 2010-03-17 - Chin)
mkdir /hive/data/genomes/felCatV17e/bed/cpgIsland
cd /hive/data/genomes/felCatV17e/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 ../../felCatV17e.2bit:$C stdout \
| maskOutFa stdin hard hardMaskedFa/${C}.fa
done
ssh swarm
cd /hive/data/genomes/felCatV17e/bed/cpgIsland
mkdir results
cut -f1 ../../chrom.sizes > chr.list
cat << '_EOF_' > template
#LOOP
./runOne $(path1) {check out exists results/$(path1).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
# in jobList:
# .....
# ./runOne chrD3 {check out exists results/chrD3.cpg}
# wc -l jobList
# .....
# 104054 jobList
# Completed: 104054 of 104054 jobs
# CPU time in finished jobs: 214s 3.57m 0.06h 0.00d 0.000 y
# IO & Wait Time: 935922s 15598.69m 259.98h 10.83d 0.030 y
# Average job time: 9s 0.15m 0.00h 0.00d
# Longest finished job: 174s 2.90m 0.05h 0.00d
# Submission to last job: 1136s 18.93m 0.32h 0.01d
# 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/felCatV17e/bed/cpgIsland
hgLoadBed felCatV17e cpgIslandExt -tab \
-sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
# Reading cpgIsland.bed
# Loaded 57950 elements of size 10
# Sorted
# Creating table definition for cpgIslandExt
# Saving bed.tab
# Loading felCatV17e
# check it with featureBits felCatV17e cpgIslandExt
# 36646772 bases of 1990635005 (1.841%) in intersection
# cleanup
rm -fr hardMaskedFa
#######################################################################
# felCatV17e Cat BLASTZ/CHAIN/NET (DONE - 2010-03-22 - Chin)
screen # use a screen to manage this multi-day job
mkdir /hive/data/genomes/canFam2/bed/lastzFelCatV17e.2010-03-22
cd /hive/data/genomes/canFam2/bed/lastzFelCatV17e.2010-03-22
cat << '_EOF_' > DEF
# dog vs. cat
# maximum M allowed with lastz is only 254
BLASTZ_M=254
# TARGET: Dog canFan2
SEQ1_DIR=/scratch/data/canFam2/canFam2.2bit
SEQ1_LEN=/scratch/data/canFam2/chrom.sizes
SEQ1_CHUNK=20000000
SEQ1_LAP=10000
SEQ1_LIMIT=5
# QUERY: Cat (felCatV17e)
SEQ2_DIR=/scratch/data/felCatV17e/felCatV17e.2bit
SEQ2_LEN=/scratch/data/felCatV17e/chrom.sizes
SEQ2_LIMIT=50
SEQ2_CHUNK=10000000
SEQ2_LAP=0
BASE=/hive/data/genomes/canFam2/bed/lastzFelCatV17e.2010-03-22
TMPDIR=/scratch/tmp
'_EOF_'
# << this line keeps emacs coloring happy
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
`pwd`/DEF \
-syntenicNet -noDbNameCheck \
-chainMinScore=3000 -chainLinearGap=medium \
-workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
> do.log 2>&1 &
# real 287m24.258s
cat fb.canFam2.chainFelCatV17eLink.txt
# 1481040604 bases of 2384996543 (62.098%) in intersection
mkdir /hive/data/genomes/felCatV17e/bed/blastz.canFam2.swap
cd /hive/data/genomes/felCatV17e/bed/blastz.canFam2.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/canFam2/bed/lastzFelCatV17e.2010-03-22/DEF \
-swap -syntenicNet -noDbNameCheck \
-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 251m14.341s
# *** All done ! Elapsed time: 251m14s
# *** Make sure that goldenPath/felCatV17e/vsCanFam2/README.txt is
# accurate.
# *** Add {chain,net}CanFam2 tracks to trackDb.ra if necessary.
cat fb.felCatV17e.chainCanFam2Link.txt
# 1467506008 bases of 1990635005 (73.720%) in intersection
#####################################################################
# LASTZ Mouse Swap (DONE - 2010-03-22 - Chin)
cd /hive/data/genomes/mm9/bed/lastzFelCatV17e.2010-03-22
cat fb.mm9.chainFelCatV17eLink.txt
# 637007193 bases of 2620346127 (24.310%) in intersection
# make sure the link exists:
# cd /hive/data/genomes/mm9/bed
# ln -s /hive/data/genomes/mm9/bed/lastzFelCatV17e.2010-03-22 \
# /lastz.felCatV17e
# swap
mkdir /hive/data/genomes/felCatV17e/bed/blastz.mm9.swap
cd /hive/data/genomes/felCatV17e/bed/blastz.mm9.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/mm9/bed/lastzFelCatV17e.2010-03-22/DEF \
-swap -syntenicNet -noDbNameCheck \
-workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 125m37.926s
cat fb.felCatV17e.chainMm9Link.txt
# 616529959 bases of 1990635005 (30.972%) in intersection
#####################################################################
# LASTZ Human Swap (DONE - 2010-03-22 - Chin)
cd /hive/data/genomes/hg19/bed/lastzFelCatV17e.2010-03-22
cat fb.hg19.chainFelCatV17eLink.txt
# 1266003011 bases of 2897316137 (43.696%) in intersection
# Swap re-run 04-07
mkdir /hive/data/genomes/felCatV17e/bed/blastz.hg19.swap
cd /hive/data/genomes/felCatV17e/bed/blastz.hg19.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/hg19/bed/lastzFelCatV17e.2010-03-22/DEF \
-swap -syntenicNet -noDbNameCheck \
-workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 200m43.245s
cat fb.felCatV17e.chainHg19Link.txt
# 1211702270 bases of 1990635005 (60.870%) in intersection
#####################################################################
# LASTZ Panda Swap (DONE - 2010-03-22 - Chin)
cd /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
cat fb.felCatV17e.chainAilMel1Link.txt
# 1503647735 bases of 1990635005 (75.536%) in intersection
# link it
# cd /hive/data/genomes/ailMel1/bed
# ln -s /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
# lastz.felCatV17e
# run rbest 2010-04-30
time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \
felCatV17e ailMel1 > rbest.log 2>&1 &
# real 234m49.556s
# re-run swap 04-26
mkdir /hive/data/genomes/felCatV17e/bed/blastz.ailMel1.swap
cd /hive/data/genomes/felCatV17e/bed/blastz.ailMel1.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22/DEF \
-swap -syntenicNet -noDbNameCheck \
-workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
# real 155m57.566s
cd /hive/data/genomes/ailMel1/bed/blastz.felCatV17e.swap
cat fb.ailMel1.chainFelCatV17eLink.txt
# 1507273252 bases of 2245312831 (67.130%) in intersection
# Due to the SEQ1 and SEQ2 spec'ed in
# /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22/DEF
# need to copy the dDirectory from ailMel1 to felCatv17e
# without this, 6way step will have problem
# cd /hive/data/genomes/ailMel1/bed/blastz.felCatV17e.swap]
# cp -pr * /hive/data/genomes/felCatV17e/bed/blastz.ailMel1.swap/
# cd /hive/data/genomes/felCatV17e/bed
# ln -s lastzAilMel1.2010-03-22 lastz.ailMel1
# mkdir /hive/data/genomes/felCatV17e/bed/lastzAilMel1.2010-03-22
# cd /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
# cp -rp * /hive/data/genomes/felCatV17e/bed/lastzAilMel1.2010-03-22/
# cd /hive/data/genomes/felCatV17e/bed
# ln -s /hive/data/genomes/felCatV17e/bed/lastzAilMel1.2010-03-22 lastz.ailMel1
# cd /hive/data/genomes/ailMel1/bed
# ln -s blastz.felCatV17e.swap lastz.felCatV17e
#####################################################################
# LASTZ Oppsom Swap (DONE - 2010-03-22 - Chin)
cd /hive/data/genomes/monDom5/bed/lastzFelCatV17e.2010-03-22
cat fb.monDom5.chainFelCatV17eLink.txt
# 178616721 bases of 3501660299 (5.101%) in intersection
# make sure the link exist
# cd /hive/data/genomes/monDom5/bed
# ln -s lastzFelCatV17e.2010-03-22 lastz.felCatV17e
# Swap
mkdir /hive/data/genomes/felCatV17e/bed/blastz.monDom5.swap
cd /hive/data/genomes/felCatV17e/bed/blastz.monDom5.swap
time nice -n +19 doBlastzChainNet.pl -verbose=2 \
/hive/data/genomes/monDom5/bed/lastzFelCatV17e.2010-03-22/DEF \
-swap -syntenicNet -noDbNameCheck \
-workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
-chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 &
# *** All done ! Elapsed time: 100m29s
# *** Make sure that goldenPath/felCatV17e/vsMonDom5/README.txt is
# accurate.
# *** Add {chain,net}MonDom5 tracks to trackDb.ra if necessary.
cat fb.felCatV17e.chainMonDom5Link.txt
# 166499264 bases of 1990635005 (8.364%) in intersection
############################################################################
# ctgPos2 track - showing clone sequence locations on chromosomes
# (DONE - 2010-03-30 - Chin)
mkdir /hive/data/genomes/felCatV17e/bed/ctgPos2
cd /hive/data/genomes/felCatV17e/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 = <FH>) {
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 775 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 felCatV17e ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql \
ctgPos2.tab
# add the track ctgPos2 to src/hg/makeDb/trackDb/cat/felCatV17e/trackDb.ra
# at src/makeDb/trackdb do "make update" or/and "make alpha"
#########################################################################
# all.joiner update, downloads and in pushQ - (DONE - 2010-04-09 Chin)
cd $HOME/kent/src/hg/makeDb/schema
# fixup all.joiner until this is a clean output
joinerCheck -database=felCatV17e -all all.joiner
mkdir /hive/data/genomes/felCatV17e/goldenPath
cd /hive/data/genomes/felCatV17e/goldenPath
makeDownloads.pl felCatV17e > do.log 2>&1
# vi /hive/data/genomes/felCatV17e/goldenPath/database/README.txt
# vi /hive/data/genomes/felCatV17e/goldenPath/bigZips/README.txt
# now ready for pushQ entry
mkdir /hive/data/genomes/felCatV17e/pushQ
cd /hive/data/genomes/felCatV17e/pushQ
makePushQSql.pl felCatV17e > felCatV17e.pushQ.sql 2> stderr.out
# check for errors in stderr.out, some are OK, e.g.:
# WARNING: felCatV17e does not have seq
# WARNING: felCatV17e does not have extFile
# WARNING: felCatV17e does not have tableDescriptions
# copy it to hgwbeta
# use hgcat in .hg.conf.beta
scp -p felCatV17e.pushQ.sql hgwbeta:/tmp
ssh hgwbeta
cd /tmp
hgsql qapushq < felCatV17e.pushQ.sql
# in that pushQ entry walk through each entry and see if the
# sizes will set properly
# reset to hguser in .hg.conf.beta
#####################################################################
-## 6-Way Multiz (working - 2010-04-19 - Chin)
+## 6-Way Multiz (DONE - 2010-04-19 - Chin)
# use /cluster/home/chinhli/kent/src/hg/utils/phyloTrees/49way.nh
mkdir /hive/data/genomes/felCatV17e/bed/multiz6way
cd /hive/data/genomes/felCatV17e/bed/multiz6way
/cluster/bin/phast/tree_doctor \
--prune-all-but=felCat3,hg19,mm9,canFam2,monDom5,ailMel1 \
--rename="felCat3 -> felCatV17e " \
/cluster/home/chinhli/kent/src/hg/utils/phyloTrees/49way.nh > 6way.nh
# rearrange felCatV17e to the top, get some help from tree_doctor:
/cluster/bin/phast/tree_doctor --name-ancestors --reroot felCatV17e \
--with-branch 6way.nh
# edit out the ancestors, and move felCatV17e from the bottom to
# the top, resulting in this tree:
(((felCatV17e:0.098612,(canFam2:0.052458,ailMel1:0.050000):0.050000):0.093470,(hg19:0.144018,mm9:0.356483):0.020593):0.258392,monDom5:0.340786);
# more rearranging after seeing what the distance table looks like
# below to get them appearing as much as possible in their
# distance order top to bottom: (TO_DO)
# use the phyloGif to get the tree image
/cluster/bin/x86_64/phyloGif -phyloGif_tree=6way.nh \
-phyloGif_width=260 -phyloGif_height=260 > 6way.gif
# more rearranging after seeing what the distance table looks like
# below to get them appearing as much as possible in their
# distance order top to bottom: (TO_DO)
# usee this specification in the phyloGif tool after changing the names:
/cluster/bin/phast/tree_doctor \
--rename="felCatV17e -> Cat ; hg19 -> Human ; ailMel1 -> Panda ; canFam2 -> Dog ; mm9 -> Mouse ; monDom5 -> Opossum " 6way.nh
# http://genome.ucsc.edu/cgi-bin/phyloGif
# to obtain a gif image for htdocs/images/phylo/felCatV17e_6way.gif
/cluster/bin/x86_64/phyloGif -phyloGif_tree=6way.nh \
-phyloGif_width=260 -phyloGif_height=260 > felCatV17e_6way.gif
/cluster/bin/phast/all_dists 6way.nh > 6way.distances.txt
# make sure all symlinks lastz.DB -> lastzDb-date
# exist here and at the swap locations, the perl script expects this
# in order to find featureBits numbers.
cd /hive/data/genomes/felCatV17e/bed]
ln -s blastz.hg19.swap lastz.hg19
ln -s blastz.canFam2.swap lastz.canFam2
ln -s blastz.mm9.swap lastz.mm9
ln -s blastz.monDom5.swap lastz.monDom5
ln -s /hive/data/genomes/felCatV17e/bed/lastzAilMel1.2010-03-22 /
lastz.ailMel1
cd /hive/data/genomes/felCatV17e/bed/multiz6way
# Use 6way.distances.txt to create the table below
# with this perl script:
cat << '_EOF_' > sizeStats.pl
#!/usr/bin/env perl
use strict;
use warnings;
open (FH, "grep -y felCatV17e 6way.distances.txt | sort -k3,3n|") or
die "can not read 6way.distances.txt";
my $count = 0;
while (my $line = <FH>) {
chomp $line;
my ($felCatV17e, $D, $dist) = split('\s+', $line);
my $chain = "chain" . ucfirst($D);
my $B="/hive/data/genomes/felCatV17e/bed/lastz.$D/fb.felCatV17e." .
$chain . "Link.txt";
my $chainLinkMeasure =
`awk '{print \$5}' ${B} 2> /dev/null | sed -e "s/(//; s/)//"`;
chomp $chainLinkMeasure;
$chainLinkMeasure = 0.0 if (length($chainLinkMeasure) < 1);
$chainLinkMeasure =~ s/\%//;
my $swapFile="/hive/data/genomes/${D}/bed/lastz.felCatV17e/fb.${D}.chainFelCatV17eLink.txt";
my $swapMeasure = "N/A";
if ( -s $swapFile ) {
$swapMeasure =
`awk '{print \$5}' ${swapFile} 2> /dev/null | sed -e "s/(//; s/)//"`;
chomp $swapMeasure;
$swapMeasure = 0.0 if (length($swapMeasure) < 1);
$swapMeasure =~ s/\%//;
}
my $orgName=
`hgsql -N -e "select organism from dbDb where name='$D';" hgcentraltest`;
chomp $orgName;
if (length($orgName) < 1) {
$orgName="N/A";
}
++$count;
++$count;
if ($swapMeasure eq "N/A") {
printf "# %02d %.4f - %s %s\t(%% %.3f) (%s)\n", $count, $dist,
$orgName, $D, $chainLinkMeasure, $swapMeasure
} else {
printf "# %02d %.4f - %s %s\t(%% %.3f) (%% %.3f)\n", $count, $dist,
$orgName, $D, $chainLinkMeasure, $swapMeasure
}
}
close (FH);
'_EOF_'
# << happy emacs
chmod +x ./sizeStats.pl
./sizeStats.pl
#
# If you can fill in all the numbers in this table, you are ready for
# the multiple alignment procedure
#
# featureBits chainLink measures
# chainFlCatV17eLink
# distance on felCatV17e on other
# 02 0.1986 - Panda ailMel1 (% 75.536) (% 67.130)
# 04 0.2011 - Dog canFam2 (% 73.720) (% 62.098)
# 06 0.3567 - Human hg19 (% 60.870) (% 43.696)
# 08 0.5692 - Mouse mm9 (% 30.972) (% 24.310)
# 10 0.7913 - Opossum monDom5 (% 8.364) (% 5.101)
# create species list and stripped down tree for autoMZ
sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \
6way.nh > tmp.nh
echo `cat tmp.nh` > tree-commas.nh
echo `cat tree-commas.nh` | sed 's/ //g; s/,/ /g' > tree.nh
sed 's/[()]//g; s/,/ /g' tree.nh > species.list
# collect the single whole mafs into one place for splitting:
mkdir singleMaf
cd singleMafs
ln -s ../../lastz.hg19/axtChain/felCatV17e.hg19.synNet.maf.gz .
ln -s ../../lastz.mm9/axtChain/felCatV17e.mm9.synNet.maf.gz .
ln -s ../../lastz.canFam2/axtChain/felCatV17e.canFam2.synNet.maf.gz .
ln -s ../../lastz.monDom5/axtChain/felCatV17e.monDom5.synNet.maf.gz .
# N50 for ailMel1 is 1281781 (less than 10 ~ 20 MB) use rbest
ln -s ../../lastz.ailMel1/mafRBestNet/felCatV17e.ailMel1.rbest.maf.gz .
# to use rbest or net, use n50.pl chrom.size to tell
mkdir /hive/data/genomes/felCatV17e/bed/multiz6way/splitMaf
cd /hive/data/genomes/felCatV17e/bed/multiz6way/splitMaf
for D in ailMel1
do
mkdir ${D}
mafSplit -useHashedName=8 -byTarget /dev/null ${D}/ \
../singleMafs/felCatV17e.${D}.rbest.maf.gz
done
for D in hg19 mm9 canFam2 monDom5
do
mkdir ${D}
mafSplit -useHashedName=8 -byTarget /dev/null ${D}/ \
../singleMafs/felCatV17e.${D}.synNet.maf.gz
done
cd /hive/data/genomes/felCatV17e/bed/multiz6way
mkdir penn
cp -p /cluster/bin/penn/multiz.2008-11-25/multiz penn
cp -p /cluster/bin/penn/multiz.2008-11-25/maf_project penn
cp -p /cluster/bin/penn/multiz.2008-11-25/maf_project penn
cp -p /cluster/bin/penn/multiz.2008-11-25/autoMZ penn
ssh swarm
mkdir /hive/data/genomes/felCatV17e/bed/multiz6way/run
mkdir /hive/data/genomes/felCatV17e/bed/multiz6way/run/maf
cd /hive/data/genomes/felCatV17e/bed/multiz6way/run
# set the db and pairs directories here
cat > autoMultiz.csh << '_EOF_'
#!/bin/csh -ef
set db = felCatV17e
set topDir = /hive/data/genomes/$db/bed/multiz6way
set c = $1
set result = $2
set pennBin = $topDir/penn
set run = `/bin/pwd`
set tmp = /scratch/tmp/$db/multiz.$c
set pairs = $topDir/splitMaf
/bin/rm -fr $tmp
/bin/mkdir -p $tmp
/bin/cp -p $topDir/tree.nh $topDir/species.list $tmp
pushd $tmp > /dev/null
foreach s (`/bin/sed -e "s/^$db //" species.list`)
set in = $pairs/$s/$c.maf
set out = $db.$s.sing.maf
if (-e $in.gz) then
/bin/zcat $in.gz > $out
if (! -s $out) then
echo "##maf version=1 scoring=autoMZ" > $out
endif
else if (-e $in) then
/bin/ln -s $in $out
else
echo "##maf version=1 scoring=autoMZ" > $out
endif
end
set path = ($pennBin $path); rehash
$pennBin/autoMZ + T=$tmp E=$db "`cat tree.nh`" $db.*.sing.maf $c.maf \
> /dev/null
popd > /dev/null
/bin/rm -f $result
/bin/cp -p $tmp/$c.maf $result
/bin/rm -fr $tmp
/bin/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+ /hive/data/genomes/felCatV17e/bed/multiz6way/run/maf/$(root1).maf}
#ENDLOOP
'_EOF_'
# << happy emacs
find ../splitMaf -type f | grep "/[0-9][0-9][0-9].maf" \
| xargs -L 1 basename | sort -u > chr.part.list
gensub2 chr.part.list single template jobList
para -ram=8g create jobList
# put the split mafs back together into a single result
head -q -n 1 maf/000.maf > felCatV17e.6way.maf
for F in maf/*.maf
do
grep -h -v "^#" ${F} >> felCatV17e.6way.maf
done
tail -q -n 1 maf/000.maf >> felCatV17e.6way.maf
- tail -q -n 1 maf/hg19_${C}.*.maf | sort -u >> ../maf/${C}.maf
# real 13m32.340s
# load tables for a look
mkdir -p /gbdb/felCatV17e/multiz6way/maf
- cd /hive/data/genomes/felCatV17e/bed/multiz6way/maf
+ # cd /hive/data/genomes/felCatV17e/bed/multiz6way/maf
+ cd /hive/data/genomes/felCatV17e/bed/multiz6way/run
+
ln -s `pwd`/felCatV17e.6way.maf \
/gbdb/felCatV17e/multiz6way/maf/multiz6way.maf
# this generates an immense multiz6way.tab file in the directory
# where it is running. Best to run this over in scratch.
cd /data/tmp
time nice -n +19 hgLoadMaf \
-pathPrefix=/gbdb/felCatV17e/multiz6way/maf felCatV17e multiz6way
- # Loaded 13316945 mafs in 1 files from
- # /gbdb/felCatV17e/multiz6way/maf
- # real 9m9.365s
+ # Indexing and tabulating /gbdb/felCatV17e/multiz6way/maf/multiz6way.maf
+ # Loading multiz6way into database
+ # Loaded 5500990 mafs in 1 files from /gbdb/felCatV17e/multiz6way/maf
+ #
+ # real 5m56.770s
# load summary table
time nice -n +19 cat /gbdb/felCatV17e/multiz6way/maf/*.maf \
| hgLoadMafSummary felCatV17e -minSize=30000 -verbose=2 \
-mergeGap=1500 -maxSize=200000 multiz6waySummary stdin
-# Created 2330531 summary blocks from 99659162 components and
-# 13316945 mafs from stdin
- # real 17m54.685s
+ # Created 1029626 summary blocks from 15177161 components and
+ # 5500990 mafs from stdin
+ # real 9m46.074s
+
+ # Gap Annotation
+ # prepare bed files with gap info
+ mkdir /hive/data/genomes/felCatV17e/bed/multiz6way/anno
+ cd /hive/data/genomes/felCatV17e/bed/multiz6way/anno
+ mkdir maf run
+
+ # most of these will already exist from previous multiple
+ # alignments
+ # remove the echo from in front of the twoBitInfo command to get
+ # them
+ # to run if this loop appears to be correct
+ for DB in `cat ../species.list`
+do
+ CDIR="/hive/data/genomes/${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 `sed -e "s/felCatV17e //" ../../species.list`
+do
+ echo "${DB} "
+ ln -s /hive/data/genomes/${DB}/${DB}.N.bed ${DB}.bed
+ echo ${DB}.bed >> nBeds
+ ln -s /hive/data/genomes/${DB}/chrom.sizes ${DB}.len
+ echo ${DB}.len >> sizes
+done
+ # the annotation step requires large memory, run on memk nodes
+
+ ssh memk
+ cd /hive/data/genomes/felCatV17e/bed/multiz6way/anno/run
+ ls ../../run/maf | sed -e "s/.maf//" > chr.list
+ cat << '_EOF_' > template
+#LOOP
+./anno.csh $(root1) {check out line+ ../maf/$(root1).maf}
+#ENDLOOP
+'_EOF_'
+ # << happy emacs
+
+ cat << '_EOF_' > anno.csh
+#!/bin/csh -fe
+
+set inMaf = ../../run/maf/$1.maf
+set outMaf = ../maf/$1.maf
+rm -f $outMaf
+mafAddIRows -nBeds=nBeds $inMaf /hive/data/genomes/felCatV17e/felCatV17e.2bit $outMaf
+'_EOF_'
+ # << happy emacs
+ chmod +x anno.csh
+
+ gensub2 chr.list single template jobList
+ para -ram=30g create jobList
+ # specify lots of ram to get one job per node
+ para -ram=30g push
+ #
+# para time
+# Completed: 256 of 256 jobs
+# CPU time in finished jobs: 3735s 62.25m 1.04h 0.04d 0.000 y
+# IO & Wait Time: 2382s 39.70m 0.66h 0.03d 0.000 y
+# Average job time: 24s 0.40m 0.01h 0.00d
+# Longest finished job: 600s 10.00m 0.17h 0.01d
+# Submission to last job: 787s 13.12m 0.22h 0.01d
+
+ ssh hgwdev
+ rm -fr /gbdb/felCatV17e/multiz6way/maf
+ mkdir /gbdb/felCatV17e/multiz6way/maf
+ cd /hive/data/genomes/felCatV17e/bed/multiz6way/anno/maf
+ ln -s `pwd`/*.maf /gbdb/felCatV17e/multiz6way/maf/
+ # by loading this into the table multiz6way, it will replace the
+ # previously loaded table with the unannotated mafs
+ # huge temp files are made, do them on local disk
+ cd /data/tmp
+ time nice -n +19 hgLoadMaf \
+ -pathPrefix=/gbdb/felCatV17e/multiz6way/maf felCatV17e multiz6way
+ # real 93m33.812s
+ # Loading multiz6way into database
+ # Loaded 6501831 mafs in 256 files from /gbdb/felCatV17e/multiz6way/maf
+
+ time nice -n +19 cat /gbdb/felCatV17e/multiz6way/maf/*.maf \
+ | hgLoadMafSummary felCatV17e -minSize=30000 -mergeGap=1500 \
+ -maxSize=200000 multiz6waySummary stdin
+ # Created 1029626 summary blocks from 15177161 components \
+ # and 6501831 mafs from stdin
+ # Loading into felCatV17e table multiz6waySummary...
+ # Loading complete
+ # real 72m12.808s
+
+ # by loading this into the table multiz6waySummary, it will
+ # replace
+ # the previously loaded table with the unannotated mafs
+ # remove the multiz6way*.tab files in this /data/tmp directory
+# -rw-rw-r-- 1 338738590 May 6 15:49 multiz6way.tab
+# -rw-rw-r-- 1 49960609 May 6 17:04 multiz6waySummary.tab
+ wc -l multiz6way*.tab
+ # 6501831 multiz6way.tab
+ # 1029626 multiz6waySummary.tab
+ # 7531457 total
+ rm multiz6way*.tab
+
+ # create some downloads
+ mkdir -p /hive/data/genomes/felCatV17e/bed/multiz6way/download/maf
+ cd /hive/data/genomes/felCatV17e/bed/multiz6way/download/maf
+ # time cp -p ../../anno/maf/chr*.maf .
+ time cp -p ../../anno/maf/*.maf .
+ # real 6m55.315s
+ time gzip --rsyncable *.maf
+ # real 30m37.474s
+ time md5sum *.gz > md5sum.txt
+ # real 1m36.814s
+ # user 0m13.210s
+ # sys 0m2.212s
####################################################################
# HUMAN (hg18) PROTEINS TRACK (DONE braney 2010-05-05)
# bash if not using bash shell already
cd /cluster/data/felCatV17e
mkdir /cluster/data/felCatV17e/blastDb
awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
twoBitToFa -seqList=1meg.lst felCatV17e.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 felCatV17e.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
# 2765
mkdir -p /cluster/data/felCatV17e/bed/tblastn.hg18KG
cd /cluster/data/felCatV17e/bed/tblastn.hg18KG
echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst
wc -l query.lst
# 2765 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/2765) = 290.143300
mkdir -p kgfa
split -l 290 /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
# 127 127 1651 kg.lst
mkdir -p blastOut
for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done
tcsh
cd /cluster/data/felCatV17e/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/felCatV17e/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/felCatV17e/bed/tblastn.hg18KG
gensub2 query.lst kg.lst blastGsub blastSpec
para create blastSpec
# para try, check, push, check etc.
para time
# Completed: 351155 of 351155 jobs
# CPU time in finished jobs: 15303249s 255054.15m 4250.90h 177.12d 0.485 y
# IO & Wait Time: 2843300s 47388.34m 789.81h 32.91d 0.090 y
# Average job time: 52s 0.86m 0.01h 0.00d
# Longest finished job: 134s 2.23m 0.04h 0.00d
# Submission to last job: 19084s 318.07m 5.30h 0.22d
ssh swarm
cd /cluster/data/felCatV17e/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: 127 of 127 jobs
# CPU time in finished jobs: 544804s 9080.07m 151.33h 6.31d 0.017 y
# IO & Wait Time: 82233s 1370.55m 22.84h 0.95d 0.003 y
# Average job time: 4937s 82.29m 1.37h 0.06d
# Longest finished job: 24063s 401.05m 6.68h 0.28d
# Submission to last job: 24074s 401.23m 6.69h 0.28d
cd /cluster/data/felCatV17e/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: 50782 failed: 0 errors: 0
# load table
ssh hgwdev
cd /cluster/data/felCatV17e/bed/tblastn.hg18KG
hgLoadPsl felCatV17e blastHg18KG.psl
# check coverage
featureBits felCatV17e blastHg18KG
# 23826621 bases of 1990635005 (1.197%) in intersection
featureBits felCatV17e blastHg18KG refGene -enrichment
# blastHg18KG 1.197%, refGene 0.021%, both 0.013%, cover 1.12%, enrich 52.73x
featureBits felCatV17e blastHg18KG xenoRefGene -enrichment
# blastHg18KG 1.197%, xenoRefGene 2.066%, both 1.010%, cover 84.36%, enrich 40.83x
rm -rf blastOut
#end tblastn