32b87a2e38fd80d8c3f0c1ffcf536f903bbf29c6 jnavarr5 Tue May 5 15:28:59 2026 -0700 Removing lines from the makedoc since these tracks are now done, and not currently running, refs #28259 diff --git src/hg/makeDb/doc/hg38/cactus447.txt src/hg/makeDb/doc/hg38/cactus447.txt index 704e2d2f32d..6520664f7fb 100644 --- src/hg/makeDb/doc/hg38/cactus447.txt +++ src/hg/makeDb/doc/hg38/cactus447.txt @@ -1,1463 +1,1462 @@ ################################################################################ #### 2023-06-13 - fetch maf.gz - Hiram mkdir -p /hive/data/genomes/hg38/bed/cactus447/singleCover cd /hive/data/genomes/hg38/bed/cactus447/singleCover obtained maf file from Glenn, from the 'courtyard' login: /nanopore/cgl/glenn-scratch/new-maf/447-mammalian-2022v1-single-copy.maf.gz -rw-r--r-- 1 hickey cgl 881994218891 Jun 5 11:06 /nanopore/cgl/glenn-scratch/new-maf/447-mammalian-2022v1-single-copy.maf.gz There was some go around with a different version. This one is 'single coverage' -rw-r--r-- 1 881994218891 Jun 5 11:06 447-mammalian-2022v1-single-copy.maf.gz ################################################################################ #### split this file into individual chromosome files 2023-06-21 - Hiram mkdir /hive/data/genomes/hg38/bed/cactus447/splitSingleCover cd /hive/data/genomes/hg38/bed/cactus447/splitSingleCover time mafSplit -byTarget -useFullSequenceName /dev/null ./ ../singleCover/447-mammalian-2022v1-single-copy.maf.gz Splitting 1 files by target sequence -- ignoring first argument /dev/null real 804m35.124s user 1061m9.580s sys 101m19.894s ################################################################################ ### convert the sequence names to hg38 chrom names 2023-06-12 to 06-14 - Hiram mkdir /hive/data/genomes/hg38/bed/cactus447/ucscNames cd /hive/data/genomes/hg38/bed/cactus447/ucscNames runOne script for processing each file: #!/bin/bash set -beEu -o pipefail export cactusMaf="${1}" export srcMaf="../splitSingleCover/${cactusMaf}" cactusId=${cactusMaf%.maf} export sedCommand=`grep -w "${cactusId}" cactus.to.ucsc.sed` export ucscName=`grep -w "${cactusId}" 2022v1.hg38.chrom.equiv.txt | cut -f2` case "${ucscName}" in chrM) if [ ${srcMaf} -nt ${ucscName}.maf ]; then sed -f cactus.to.ucsc.sed ${srcMaf} > ${ucscName}.maf touch -r ${srcMaf} ${ucscName}.maf fi ;; chr[0-9][0-9]) rm -f "${ucscName}.maf" ln ${srcMaf} "${ucscName}.maf" ;; chr[0-9]) rm -f "${ucscName}.maf" ln ${srcMaf} "${ucscName}.maf" ;; G*.maf) K*.maf) if [ ${srcMaf} -nt ${ucscName}.maf ]; then sed -e "${sedCommand}" ${srcMaf} > ${ucscName}.maf touch -r ${srcMaf} ${ucscName}.maf else printf "DONE: ${ucscName}.maf\n" 1>&2 fi ;; esac if [ "${ucscName}.maf" -nt "nameScan/${ucscName}.seqCount.txt" ]; then grep "^s." "${ucscName}.maf" | cut -d' ' -f2 | sort | uniq -c | sort -rn \ > "nameScan/${ucscName}.seqCount.txt" touch -r "${ucscName}.maf" "nameScan/${ucscName}.seqCount.txt" sed -e "s/^ \+//;" "nameScan/${ucscName}.seqCount.txt" | cut -d' ' -f2 \ | cut -d'.' -f1 | sort | uniq -c \ | sort -rn > "nameScan/${ucscName}.species.txt" fi ##################################################### template: #LOOP runOne $(path1) #ENDLOOP using list of maf files: ls ../splitSingleCover | grep maf > cactus.maf.list gensub2 cactus.maf.list single template jobList ### made mistake of running too many of those at the same time ### and their use of large amounts of memory caused hgwdev to ### hang up and needed a reboot to recover. So, run these carefully ### only a couple at a time on hgwdev and on hgcompute-01, the large ### ones get as large as 800 Gb in memory. ### get the species names from the nameScan: sed -e 's/^ \+//;' nameScan/chr2.species.txt | cut -d' ' -f2 \ | sort > ../species.list.txt ################################################################################ ### get the chrom.sizes out of the maf files 2023-06-19 - Hiram mkdir /hive/data/genomes/hg38/bed/cactus447/chromSizes cd /hive/data/genomes/hg38/bed/cactus447/chromSizes ### runOne script to run each one: #!/bin/bash export inMaf="${1}" export result="${2}" grep "^s " "${inMaf}" | awk '{printf "%s\t%s\n", $2,$6}' | sort -u > "${result} ### working on the maf list: ls ../ucscNames/chr*.maf > maf.list ### template to construct jobList #LOOP runOne $(path1) {check out line+ result/$(root1).txt} #ENDLOOP gensub2 maf.list single template jobList para -ram=6g create jobList Completed: 101 of 101 jobs CPU time in finished jobs: 90801s 1513.35m 25.22h 1.05d 0.003 y IO & Wait Time: 0s 0.00m 0.00h 0.00d 0.000 y Average job time: 762s 12.69m 0.21h 0.01d Longest finished job: 7952s 132.53m 2.21h 0.09d Submission to last job: 8022s 133.70m 2.23h 0.09d ### sizes.pl script to get the answer out of the result/*.txt files: #!/usr/bin/env perl use strict; use warnings; my %chromSizes; # key is assembly name, value is a pointer to a hash # with key chromName, value size open (my $FH, "-|", "ls result/*.txt") or die "can not read ls result/*.txt"; while (my $file = <$FH>) { chomp $file; printf STDERR "# reading: %s\n", $file; open (my $R, "<", $file) or die "can not read $file"; while (my $line = <$R>) { chomp $line; my @a = split('\.', $line); if (! defined($chromSizes{$a[0]})) { my %h; $chromSizes{$a[0]} = \%h; } my $asm = $chromSizes{$a[0]}; if (scalar(@a) == 3) { my ($a2, $size) = split('\s+', $a[2]); my $chr = "$a[1].$a2"; if (defined($asm->{$chr})) { if ($asm->{$chr} != $size) { printf STDERR "ERROR: size mismatch: %d != %d for %s\t%s.%s\n", $asm->{$chr}, $size, $a[0], $a[1], $a[2]; exit 255; } } else { $asm->{$chr} = $size; } # printf "%s\t%s.%s\n", $a[0], $a[1], $a[2]; } elsif (scalar(@a) == 2) { my ($chr, $size) = split('\s+', $a[1]); if (defined($asm->{$chr})) { if ($asm->{$chr} != $size) { printf STDERR "ERROR: size mismatch: %d != %d for %s\t%s\n", $asm->{$chr}, $size, $a[0], $a[1]; exit 255; } } else { $asm->{$chr} = $size; } # printf "%s\t%s\n", $a[0], $a[1]; } else { printf STDERR "ERROR: not 2 or 3: '%s'\n", $line; exit 255; } } close ($R); } close ($FH); foreach my $asm (sort keys %chromSizes) { printf STDERR "# output %s\n", $asm; my $out = "sizes/${asm}.chrom.sizes"; open ($FH, "|-", "sort -k2,2nr > $out") or die "can not write to $out"; my $h = $chromSizes{$asm}; foreach my $chr (sort keys %$h) { # printf "%s\t%s\t%d\n", $asm, $chr, $h->{$chr}; printf $FH "%s\t%d\n", $chr, $h->{$chr}; } close ($FH); } mkdir sizes time (./sizes.pl) > sizes.log 2>&1 real 10m11.808s user 6m40.128s sys 3m37.475s ################################################################################ ### add the iRow annotation 2023-06-19 to 2023-06-26 - Hiram mkdir /hive/data/genomes/hg38/bed/cactus447/iRows cd /hive/data/genomes/hg38/bed/cactus447/iRows # obtained the N.bed files from Glenn ## symLink them into this directory: ls /hive/data/genomes/hg38/bed/cactus447/Nbeds/*.N.bed | grep -v "hg38.N.bed" \ | while read P do B=`basename $P | sed -e 's/N.bed/bed/;'` rm -f ./${P} ln -s "${P}" ./${B} done ### and symlink in the chrom.sizes files as .len files ls /hive/data/genomes/hg38/bed/cactus447/chromSizes/sizes/*.chrom.sizes \ | grep -v hg38.chrom.sizes | while read P do S=`basename ${P} | sed -e 's/chrom.sizes/len/;'` printf "%s\n" "${S}" ln -s "${P}" ./${S} done ln -s /hive/data/genomes/hg38/chrom.sizes hg38.len ls *.bed > nBeds ls *.len > sizes ### list of maf file to process: ls ../ucscNames/chr*.maf > maf.list ### template file to construct jobList #LOOP mafAddIRows -nBeds=nBeds $(path1) /hive/data/genomes/hg38/hg38.2bit result/$(file1) #ENDLOOP gensub2 maf.list single template jobList ### took a week to get these all done. ### convert each file to a bigMaf file: mkdir /hive/data/genomes/hg38/bed/cactus447/iRows/bigMaf cd /hive/data/genomes/hg38/bed/cactus447/iRows/bigMaf ### tried to sort them, but made a mistake on the sort -k argument: for F in ../result/chr*.maf do B=`basename $F | sed -e 's/.maf//;'` mafToBigMaf hg38 "${F}" stdout | $HOME/bin/x86_64/gnusort -k2,2 -S1000G \ --parallel=32 > "${B}.bigMaf" done real 1588m59.209s user 1259m38.364s sys 371m13.448s ### put these all together: ls -S chr*.bigMaf | while read M do cat "${M}" done > ../hg38.cactus447.bigMaf real 322m25.741s user 0m13.577s sys 82m41.509s ### since they were sorted incorrectly above, sort this one huge file: $HOME/bin/x86_64/gnusort -S500G --parallel=64 -k1,1 -k2,2n \ hg38.cactus447.bigMaf > hg38.sorted447way.bigMaf real 752m12.425s user 47m34.825s sys 371m7.394s ################################################################################ ### make the bigBed file from the sorted bigMaf 2023-06-27 bedToBigBed -verbose=2 -itemsPerSlot=4 -type=bed3+1 \ -as=$HOME/kent/src/hg/lib/bigMaf.as -tab iRows/hg38.sorted447way.bigMaf \ /hive/data/genomes/hg38/chrom.sizes hg38.cactus447way.bb # pid=61008: VmPeak: 5682928 kB # real 4157m14.091s # that time is 69 hours ############################################################################## # liftOver files to convert sequence names (DONE - 2023-08-04 - Hiram) # # The sequence names in the assemblies in the alignment do not match # the sequence names in the corresponding genome browsers assemblies mkdir /hive/data/genomes/hg38/bed/cactus447/liftOver cd /hive/data/genomes/hg38/bed/cactus447/liftOver mkdir lifts ln -s ../idKeys/447way.equivalent.txt . # using this script to construct liftOver files #!/bin/bash TOP="/hive/data/genomes/hg38/bed/cactus447/liftOver" cd "${TOP}" mkdir -p lifts cat 447way.equivalent.txt | while read L do seqName=`printf "%s" "${L}" | cut -f1` asmId=`printf "%s" "${L}" | cut -f2` seqIdKeys="/hive/data/genomes/hg38/bed/cactus447/idKeys/${seqName}/${seqName}.idKeys.txt" case ${asmId} in GCA_* | GCF_* ) gcX="${asmId:0:3}" d0="${asmId:4:3}" d1="${asmId:7:3}" d2="${asmId:10:3}" ab="/hive/data/genomes/asmHubs/allBuild/${gcX}/${d0}/${d1}/${d2}/${asmId}" buildDir=`realpath "${ab}"` asmKeys="${buildDir}/trackData/idKeys/${asmId}.idKeys.txt" asmSizes="${buildDir}/${asmId}.chrom.sizes" # asmSizes="${buildDir}/idKeys/chrom.sizes" ls -og "${asmKeys}" ;; [a-zA-Z0-9]*) asmKeys="/hive/data/genomes/${asmId}/bed/idKeys/${asmId}.idKeys.txt" asmSizes="/hive/data/genomes/${asmId}/chrom.sizes" ;; *) printf "ERROR not recognized $asmId\n" 1>&2 exit 255 ;; esac seqKeys="/hive/data/genomes/hg38/bed/cactus447/idKeys/${seqName}/${seqName}.idKeys.txt" seq2Bit="/hive/data/genomes/hg38/bed/cactus447/2bitFiles/${seqName}.2bit" if [ ! -s "${seqKeys}" ]; then printf "can not find seq idKeys\n%s\n" "${seqKeys}" 1>&2 exit 255 fi if [ ! -s "${asmKeys}" ]; then printf "can not find asm idKeys\n%s\n" "${asmKeys}" 1>&2 exit 255 fi printf "%s\t%s\n" "${asmId}" "${seqName}" 1>&2 printf "lifts/${seqName}.${asmId}.lift\n" 1>&2 printf "$seqKeys\n" 1>&2 printf "$asmKeys\n" 1>&2 printf "$asmSizes\n" 1>&2 if [ ! -s "lifts/${seqName}.${asmId}.lift" ]; then join -t$'\t' "${seqKeys}" "${asmKeys}" | sort -k3,3 \ | join -t$'\t' -2 3 <(sort -k1,1 ${asmSizes}) - \ | awk -F$'\t' '{printf "0\t%s\t%d\t%s\t%d\n", $1, $2, $4, $2}' \ > lifts/${seqName}.${asmId}.lift fi done # a couple of those ended up with duplicate names in the liftOver # due to duplicate contigs in the assemblies. They were edited # manually, one was really bad and needed this perl script to clean it: #!/usr/bin/env perl use strict; use warnings; my %idDone; # key is sequence name, value is string of names used open (my $fh, "<", "lifts/Otolemur_garnettii.otoGar3.lift") or die "can not read lifts/Otolemur_garnettii.otoGar3.lift"; while (my $line = <$fh>) { chomp $line; my @a = split('\t', $line); if (!defined($idDone{$a[1]})) { $idDone{$a[1]} = $a[3]; printf "%s\n", $line; } else { $idDone{$a[1]} .= "," . $a[3]; } } close ($fh); foreach my $seq (sort keys %idDone) { printf STDERR "%s\t%s\n", $seq, $idDone{$seq}; } ############################################################################## # MAF FRAMES (working - 2023-08-01 - Hiram) ssh hgwdev mkdir /hive/data/genomes/hg38/bed/cactus447/frames cd /hive/data/genomes/hg38/bed/cactus447/frames mkdir genes # survey all the genomes to find out what kinds of gene tracks they have printf '#/bin/bash egrep -v "GCA_|GCF_" ../idKeys/447way.equivalent.txt | cut -f2 | while read db do printf "# ${db}:" for table in `hgsql -N -e "show tables;" $db | egrep "Gene|ncbiRefSeq" | egrep -v "Curated|Cds|Link|Other|PepTable|Predicted|Psl|wgEncode|gencode|uuGene|ToGeneName|ccds|encode|Ext|Old|Pep|Select|Bak|Verify|orfeome|Multi|transMap|Pfam|knownGeneMrna|sgpGene"` do count=`hgsql -N -e "select count(*) from $table;" $db` printf " %%s: %%s" "${table}" "${count}" done printf "\\n" done ' > showGenes.sh chmod +x ./showGenes.sh ./showGenes.sh > gene.survey.txt 2>&1 & # most of them have ncbiRefSeq grep ncbiRefSeq gene.survey.txt | cut -d' ' -f2 | sed -e 's/://;' \ | sort | xargs echo | fold -s -w 72 # bosMut1 canFam4 chiLan1 colAng1 dipOrd2 eleEdw1 felCat8 fukDam1 micOch1 # mm10 myoBra1 myoDav1 myoLuc2 odoRosDiv1 orcOrc1 otoGar3 proCoq1 pteAle1 # rheMac10 sorAra2 speTri2 susScr3 tarSyr2 tupChi1 # a couple do not: grep -v ncbiRefSeq gene.survey.txt | cut -d' ' -f2 | sed -e 's/://;' \ | sort | xargs echo | fold -s -w 72 # cavApe1 eidHel1 ptePar1 # but they do have ensGene or augustusGene: egrep "cavApe1|eidHel1|ptePar1" gene.survey.txt # cavApe1: augustusGene: 10114 ensGene: 22510 xenoRefGene: 322224 # eidHel1: augustusGene: 49674 # ptePar1: augustusGene: 68613 # 1. ncbiRefSeq - add hg38 here manually since it didn't come in above # for db in hg38 # do # asmName="${db}" for db in bosMut1 canFam4 chiLan1 colAng1 dipOrd2 eleEdw1 felCat8 fukDam1 micOch1 mm10 myoBra1 myoDav1 myoLuc2 odoRosDiv1 orcOrc1 otoGar3 proCoq1 pteAle1 rheMac10 sorAra2 speTri2 susScr3 tarSyr2 tupChi1 do asmName=`grep -w "${db}" ../idKeys/447way.equivalent.txt | cut -f1` hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ncbiRefSeq" ${db} \ | genePredSingleCover stdin stdout | gzip -2c \ > genes/${asmName}.gp.gz echo -n "# ${asmName}: " genePredCheck -db=${db} genes/${asmName}.gp.gz done # hg38: checked: 22882 failed: 0 # Bos_mutus: checked: 20457 failed: 0 # CanFam4: checked: 21143 failed: 0 # Chinchilla_lanigera: checked: 20482 failed: 0 # Colobus_angolensis: checked: 20042 failed: 0 # Dipodomys_ordii: checked: 19738 failed: 0 # Elephantulus_edwardii: checked: 20384 failed: 0 # Felis_catus: checked: 20051 failed: 0 # Fukomys_damarensis: checked: 19811 failed: 0 # Microtus_ochrogaster: checked: 20048 failed: 0 # Mus_musculus: checked: 23375 failed: 0 # Myotis_brandtii: checked: 19949 failed: 0 # Myotis_davidii: checked: 18815 failed: 0 # Myotis_lucifugus: checked: 19895 failed: 0 # Odobenus_rosmarus: checked: 19346 failed: 0 # Orcinus_orca: checked: 19136 failed: 0 # Otolemur_garnettii: checked: 19536 failed: 0 # Propithecus_coquerelli: checked: 19356 failed: 0 # Pteropus_alecto: checked: 18326 failed: 0 # Macaca_mulatta: checked: 21021 failed: 0 # Sorex_araneus: checked: 19160 failed: 0 # Ictidomys_tridecemlineatus: checked: 19892 failed: 0 # Sus_scrofa: checked: 24180 failed: 0 # Carlito_syrichta: checked: 19968 failed: 0 # Tupaia_chinensis: checked: 21047 failed: 0 # 2. ensGene for db in cavApe1 do asmName=`grep -w "${db}" ../idKeys/447way.equivalent.txt | cut -f1` hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ensGene" ${db} \ | genePredSingleCover stdin stdout | gzip -2c \ > /dev/shm/${db}.tmp.gz mv /dev/shm/${db}.tmp.gz genes/${asmName}.gp.gz echo -n "# ${asmName}: " genePredCheck -db=${db} genes/${asmName}.gp.gz done # Cavia_aperea: checked: 14182 failed: 0 # 3. augustusGene for db in eidHel1 ptePar1 do asmName=`grep -w "${db}" ../idKeys/447way.equivalent.txt | cut -f1` hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from augustusGene" ${db} \ | genePredSingleCover stdin stdout | gzip -2c \ > /dev/shm/${db}.tmp.gz mv /dev/shm/${db}.tmp.gz genes/${asmName}.gp.gz echo -n "# ${asmName}: " genePredCheck -db=${db} genes/${asmName}.gp.gz done # Eidolon_helvum: checked: 31495 failed: 0 # Pteronotus_parnellii: checked: 47678 failed: 0 # verify counts for genes are reasonable: for T in genes/*.gz do printf "# $T: " zcat $T | cut -f1 | sort -u | wc -l done # genes/Bos_mutus.gp.gz: 20457 # genes/CanFam4.gp.gz: 21143 # genes/Canis_lupus_familiaris.gp.gz: 19812 # genes/Carlito_syrichta.gp.gz: 19968 # genes/Cavia_aperea.gp.gz: 14182 # genes/Cercocebus_atys.gp.gz: 20534 # genes/Chinchilla_lanigera.gp.gz: 20481 # genes/Colobus_angolensis.gp.gz: 20042 # genes/Condylura_cristata.gp.gz: 17842 # genes/Dipodomys_ordii.gp.gz: 19738 ... # genes/Pteropus_vampyrus.gp.gz: 19485 # genes/Rattus_norvegicus.gp.gz: 22854 # genes/Rhinopithecus_roxellana.gp.gz: 21027 # genes/Saimiri_boliviensis.gp.gz: 20933 # genes/Sorex_araneus.gp.gz: 19160 # genes/Sus_scrofa.gp.gz: 24043 # genes/Theropithecus_gelada.gp.gz: 20592 # genes/Tupaia_chinensis.gp.gz: 21047 # And, can pick up ncbiRefSeq for GCF equivalent assemblies for those # outside the UCSC database set of genomes grep GCF_ ../idKeys/447way.equivalent.txt | while read dbGCF do asmName=`printf "%s" "${dbGCF}" | cut -f1` asmId=`printf "%s" "${dbGCF}" | cut -f2` gcX="${asmId:0:3}" d0="${asmId:4:3}" d1="${asmId:7:3}" d2="${asmId:10:3}" buildDir="/hive/data/genomes/asmHubs/refseqBuild/$gcX/$d0/$d1/$d2/${asmId}" if [ ! -d "${buildDir}" ]; then printf "not found: %s\n" "${buildDir}" else ncbiRefSeqGp="${buildDir}/trackData/ncbiRefSeq/process/${asmId}.ncbiRefSeq.gp" sizes="${buildDir}/${asmId}.chrom.sizes" if [ -s "${ncbiRefSeqGp}" ]; then genePredSingleCover "${ncbiRefSeqGp}" stdout | gzip -c > genes/${asmName}.gp.gz printf "# %s: " "${asmName}" genePredCheck -chromSizes="${sizes}" genes/${asmName}.gp.gz else printf "MISSING: ${ncbiRefSeqGp}\n" 1>&2 fi fi done # verify counts for genes are reasonable: for T in genes/*.gz do printf "# $T: " zcat $T | cut -f1 | sort -u | wc -l done # genes/Bos_mutus.gp.gz: 20457 # genes/CanFam4.gp.gz: 21143 # genes/Canis_lupus_familiaris.gp.gz: 19812 # genes/Carlito_syrichta.gp.gz: 19968 # genes/Cavia_aperea.gp.gz: 14182 # genes/Cercocebus_atys.gp.gz: 20534 # genes/Chinchilla_lanigera.gp.gz: 20481 # genes/Colobus_angolensis.gp.gz: 20042 # genes/Condylura_cristata.gp.gz: 17842 # genes/Dipodomys_ordii.gp.gz: 19738 ... # genes/Pteropus_vampyrus.gp.gz: 19485 # genes/Rattus_norvegicus.gp.gz: 22854 # genes/Rhinopithecus_roxellana.gp.gz: 21027 # genes/Saimiri_boliviensis.gp.gz: 20933 # genes/Sorex_araneus.gp.gz: 19160 # genes/Sus_scrofa.gp.gz: 24043 # genes/Theropithecus_gelada.gp.gz: 20592 # genes/Tupaia_chinensis.gp.gz: 21047 # these gene predictions need lifting to change the sequence names mv genes genes.seqToAsm mkdir genes cd genes.seqToAsm for F in *.gp.gz do seqName=`echo $F | sed -e 's/.gp.gz//;'` liftFile="`ls -d /hive/data/genomes/hg38/bed/cactus447/liftOver/lifts/${seqName}.*.lift 2> /dev/null`" if [ -s "${liftFile}" ]; then printf "%s\n" "${seqName}" 1>&2 liftUp -type=.genepred stdout "${liftFile}" warn "${F}" | gzip -c \ > "../genes/${seqName}.gp.gz" else printf "can not find lift file for $seqName\n" 1>&2 fi done # kluster job to annotate each maf file screen -S hg38 # manage long running procedure with screen ssh ku cd /hive/data/genomes/hg38/bed/cactus447way/frames printf '#!/bin/bash set -beEu -o pipefail export C="${1}" export G="${2}" cat ../iRows/result/${C}.maf | genePredToMafFrames hg38 stdin stdout \ ${G} genes/${G}.gp.gz | gzip > parts/${C}.${G}.mafFrames.gz ' > runOne chmod +x runOne ls ../iRows/result | grep maf | sed -e "s/.maf//" > chr.list ls genes | sed -e "s/.gp.gz//" > gene.list printf '#LOOP runOne $(root1) $(root2) {check out exists+ parts/$(root1).$(root2).mafFrames.gz} #ENDLOOP ' > template mkdir parts # run on the ku kluster gensub2 chr.list gene.list template jobList para -ram=128g create jobList para try ... check ... push # Completed: 4444 of 4444 jobs # CPU time in finished jobs: 3955725s 65928.75m 1098.81h 45.78d 0.125 y # IO & Wait Time: 1735948s 28932.46m 482.21h 20.09d 0.055 y # Average job time: 1281s 21.35m 0.36h 0.01d # Longest finished job: 34637s 577.28m 9.62h 0.40d # Submission to last job: 55372s 922.87m 15.38h 0.64d # collect all results into one file: cd /hive/data/genomes/hg38/bed/cactus447way/frames find ./parts -type f | while read F; do zcat ${F} | awk '11 == NF' done | sort -k1,1 -k2,2n > cactus447wayFrames.bed bedToBigBed -type=bed3+7 -as=$HOME/kent/src/hg/lib/mafFrames.as \ -tab cactus447wayFrames.bed ../../../chrom.sizes cactus447wayFrames.bb bigBedInfo cactus447wayFrames.bb | sed -e 's/^/# /;' # version: 4 # fieldCount: 11 # hasHeaderExtension: yes # isCompressed: yes # isSwapped: 0 # extraIndexCount: 0 # itemCount: 16,063,019 # primaryDataSize: 136,782,205 # primaryIndexSize: 1,013,720 # zoomLevels: 10 # chromCount: 82 # basesCovered: 65,486,909 # meanDepth (of bases covered): 22.709888 # minDepth: 1.000000 # maxDepth: 44.000000 # std of depth: 18.692097 # -rw-rw-r-- 1 1192312943 Aug 5 15:28 cactus447wayFrames.bed # -rw-rw-r-- 1 166572857 Aug 5 15:29 cactus447wayFrames.bb gzip cactus447wayFrames.bed # verify there are frames on everything expected. Since the # frames on the HL* asemblies did not work due to sequence name # differences, this won't be everything. # (ls genes | wc shows 44 'expected') zcat cactus447wayFrames.bed.gz | awk '{print $4}' | sort | uniq -c \ | sed -e 's/^/# /;' > species.check.list wc -l species.check.list # 44 # 391443 Bos_mutus # 414721 CanFam4 # 417196 Canis_lupus_familiaris # 328151 Carlito_syrichta # 255097 Cavia_aperea # 373053 Cercocebus_atys # 398356 Chinchilla_lanigera # 350442 Colobus_angolensis # 346524 Condylura_cristata # 363333 Dipodomys_ordii ... # 379732 Pteropus_alecto # 387484 Pteropus_vampyrus # 389056 Rattus_norvegicus # 360929 Rhinopithecus_roxellana # 366019 Saimiri_boliviensis # 336296 Sorex_araneus # 402603 Sus_scrofa # 357418 Theropithecus_gelada # 372558 Tupaia_chinensis # 195495 hg38 # load the resulting file, merely for measurement purposes here ssh hgwdev cd /hive/data/genomes/hg38/bed/cactus447way/frames time hgLoadMafFrames hg38 cactus447wayFrames cactus447wayFrames.bed.gz # real 3m14.979s hgsql -e 'select count(*) from cactus447wayFrames;' hg38 # +----------+ # | count(*) | # +----------+ # | 16063019 | # +----------+ time featureBits -countGaps hg38 cactus447wayFrames # 65486909 bases of 3299210039 (1.985%) in intersection # real 1m50.539s # enable the trackDb entries: # frames https://hgdownload.soe.ucsc.edu/goldenPath/hg38/cactus447way/cactus447wayFrames.bb # irows on # zoom to base level in an exon to see codon displays # appears to work OK # do not need this loaded table: hgsql hg38 -e 'drop table cactus447wayFrames;' ######################################################################### # Try loading this track data into database to see if download files # such as upstream gene mafs can be constructed mkdir -p /gbdb/hg38/cactus447way/maf cd /hive/data/genomes/hg38/bed/cactus447/iRows/result ln -s `pwd`/chr*.maf /gbdb/hg38/cactus447way/maf/ # this generates an immense cactus447way.tab file in the directory # where it is running. Best to run this over in scratch. # This is going to take all day. cd /dev/shm time hgLoadMaf -pathPrefix=/gbdb/hg38/cactus447way/maf hg38 cactus447way # -rw-rw-r-- 1 55422976 Aug 8 09:22 cactus447way.tab hgsql -e 'select count(*) from cactus447way;' hg38 # Loading cactus447way into database # Loaded 108225530 mafs in 101 files from /gbdb/hg38/cactus447way/maf # real 440m30.144s # user 404m46.183s # sys 19m27.664s ######################################################################### # Phylogenetic tree from 447-way (DONE - 2022-08-26 - Hiram) mkdir /hive/data/genomes/hg38/bed/cactus447/4d cd /hive/data/genomes/hg38/bed/cactus447/4d # the annotated maf's are in: ../iRows/result/chr*.maf # using ncbiRefSeq for hg38, only transcribed genes and nothing # from the randoms and other misc. hgsql -Ne "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ncbiRefSeq where cdsEnd > cdsStart;" hg38 \ | egrep -E -v "chrM|chrUn|random|_alt|_fix" > ncbiRefSeq.gp wc -l *.gp # 130420 ncbiRefSeq.gp # verify it is only on the chroms: cut -f2 ncbiRefSeq.gp | sort | uniq -c | sort -rn | sed -e 's/^/ # /;' # 12841 chr1 # 9518 chr2 # 8490 chr3 # 7521 chr19 # 7519 chr11 # 7329 chr17 # 7287 chr12 # 6155 chr6 # 6151 chr10 # 5905 chr7 # 5561 chr5 # 5504 chr9 # 5351 chr4 # 5235 chr16 # 4964 chr8 # 4296 chrX # 4255 chr15 # 3959 chr14 # 3089 chr20 # 2721 chr22 # 2521 chr18 # 2471 chr13 # 1411 chr21 # 366 chrY genePredSingleCover ncbiRefSeq.gp stdout | sort > ncbiRefSeqNR.gp wc -l ncbiRefSeqNR.gp # 19542 ncbiRefSeqNR.gp ssh ku mkdir /hive/data/genomes/hg38/bed/cactus447/4d/run cd /hive/data/genomes/hg38/bed/cactus447/4d/run mkdir ../mfa # newer versions of msa_view have a slightly different operation # the sed of the gp file inserts the reference species in the chr name # these processes are going to be copying the whole chromosome maf file # to /scratch/tmp/ - these are huge files, make sure /scratch/tmp/ # on the kluster nodes is clean enough to manage the largest chromosome # maf file, for example: # -rw-rw-r-- 1 661829442743 Jun 20 18:55 chr1.maf # -rw-rw-r-- 1 682410918389 Jun 21 16:21 chr2.maf cat << '_EOF_' > 4d.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2018-03-29/bin set r = "/hive/data/genomes/hg38/bed/cactus447" set c = $1 set infile = $r/iRows/result/$2 set outfile = $3 cd /scratch/tmp # 'clean' maf, removes all chrom names, leaves only the db name perl -wpe 's/^s ([^.]+)\.\S+/s $1/' $infile > $c.maf awk -v C=$c '$2 == C {print}' $r/4d/ncbiRefSeqNR.gp | sed -e "s/\t$c\t/\thg38.$c\t/" > $c.gp set NL=`wc -l $c.gp| gawk '{print $1}'` if ("$NL" != "0") then $PHASTBIN/msa_view --4d --features $c.gp -i MAF $c.maf -o SS > $c.ss $PHASTBIN/msa_view -i SS --tuple-size 1 $c.ss > $r/4d/run/$outfile else echo "" > $r/4d/run/$outfile endif rm -f $c.gp $c.maf $c.ss '_EOF_' # << happy emacs chmod +x 4d.csh ls -1S /hive/data/genomes/hg38/bed/cactus447/iRows/result/chr*.maf \ | sed -e "s#.*cactus447/iRows/result/##" \ | egrep -E -v "chrM|chrUn|random|_alt|_fix" > maf.list printf '#LOOP 4d.csh $(root1) $(path1) {check out line+ ../mfa/$(root1).mfa} #ENDLOOP ' > template gensub2 maf.list single template jobList para -ram=128g create jobList para try ... check ... push ... etc... para time # CPU time in finished jobs: 290959s 4849.31m 80.82h 3.37d 0.009 y # IO & Wait Time: 17027s 283.79m 4.73h 0.20d 0.001 y # Average job time: 15399s 256.65m 4.28h 0.18d # Longest finished job: 36250s 604.17m 10.07h 0.42d # Submission to last job: 36265s 604.42m 10.07h 0.42d # combine mfa files ssh hgwdev cd /hive/data/genomes/hg38/bed/cactus447/4d # verify no tiny files: ls -og mfa | sort -k3nr | tail -3 # -rw-rw-r-- 1 2583636 Aug 8 10:30 chr21.mfa # -rw-rw-r-- 1 1759815 Aug 8 09:47 chrY.mfa #want comma-less species.list time /cluster/bin/phast.build/cornellCVS/phast.2018-03-29/bin/msa_view \ --aggregate "`cat ../species.list.txt`" mfa/*.mfa | sed s/"> "/">"/ \ > 4d.all.mfa # real 0m9.416s # -rw-rw-r-- 1 303080280 Aug 8 21:48 4d.all.mfa # check they are all in there: grep "^>" 4d.all.mfa | wc -l # 447 sed 's/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \ ../hg38.447way.nh sed 's/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \ ../hg38.447way.nh > tree-commas.nh # use phyloFit to create tree model (output is phyloFit.mod) time /cluster/bin/phast.build/cornellCVS/phast.2018-03-29/bin/phyloFit \ --EM --precision MED --msa-format FASTA --subst-mod REV \ --tree tree-commas.nh 4d.all.mfa # real 880m32.863s mv phyloFit.mod all.mod grep TREE all.mod # TREE: (((((((((((((hg38:0.00899231,(Pan_paniscus:0.00233221, # Pan_troglodytes:0.00214362):0.00450706):0.00200006, # (Gorilla_gorilla:0.00120794, # Gorilla_beringei:0.00131286):0.00731999):0.00848699, # (Pongo_abelii:0.00175881,Pongo_pygmaeus:0.0021602):0.015206):0.00273536, # (((((Hylobates_lar:0.000904688,Hylobates_pileatus:0.000921279):0.00271202, # (((Hylobates_abbotti:0.00182513,Hylobates_muelleri:0.00188051):0.000678509, # Hylobates_agilis:0.00217126):0.00068124, # ... # (Ceratotherium_simum:0.000693402, # Ceratotherium_simum_cottoni:0.000665708):0.00612505):0.00973125):0.0392953):0.0100285):0.0323073):0.00299066):0.00383198):0.00906124):0.0213096):0.012194, # (((Dasypus_novemcinctus:0.0705749,(Tolypeutes_matacus:0.036879, # Chaetophractus_vellerosus:0.0367115):0.0193664):0.0421356, # ((Tamandua_tetradactyla:0.0243613, # Myrmecophaga_tridactyla:0.0216596):0.0822034, # (Choloepus_didactylus:0.0057887, # Choloepus_hoffmanni:0.00621468):0.0702582):0.018824):0.0529836, # (((Trichechus_manatus:0.0621184,(Procavia_capensis:0.0107201, # Heterohyrax_brucei:0.0105857):0.149548):0.00440411, # Loxodonta_africana:0.0800445):0.0234193, # ((((Microgale_talazaci:0.118008,Echinops_telfairi:0.077048):0.158928, # Chrysochloris_asiatica:0.138211):0.0147081, # Elephantulus_edwardii:0.20827):0.00712897, # Orycteropus_afer:0.112233):0.00664842):0.0534765):0.012194); # compare these calculated lengths to what we started with /cluster/bin/phast/all_dists ../hg38.447way.nh | grep hg38 \ | sed -e "s/hg38.//;" | sort > original.dists grep TREE all.mod | sed -e 's/TREE: //;' \ | /cluster/bin/phast/all_dists /dev/stdin | grep hg38 \ | sed -e "s/hg38.//;" | sort > hg38.dists # printing out the 'original', the 'new' the 'difference' and # percent difference/delta join original.dists hg38.dists | awk '{ printf "#\t%s\t%8.6f\t%8.6f\t%8.6f\t%8.6f\n", $1, $2, $3, $2-$3, 100*($2-$3)/$3 }' | sort -k4n | cut -f2,3,4,6 | sort -k4,4n Chiropotes_albinasus 0.044221 0.153018 -71.100786 Hapalemur_gilberti 0.091144 0.238983 -61.861722 Hapalemur_meridionalis 0.091511 0.239866 -61.849116 Hapalemur_occidentalis 0.090759 0.237751 -61.826028 Hapalemur_griseus 0.091379 0.239368 -61.824889 Hapalemur_alaotrensis 0.091149 0.238685 -61.812012 Eulemur_rufus 0.090624 0.236415 -61.667407 Eulemur_fulvus 0.090508 0.236040 -61.655652 Lemur_catta 0.090748 0.236490 -61.627130 ... Tolypeutes_matacus 0.518057 0.332742 55.693300 Chaetophractus_vellerosus 0.520177 0.332575 56.408930 Catagonus_wagneri 0.578187 0.366272 57.857275 Megaderma_lyra 0.575967 0.360862 59.608659 Spermophilus_dauricus 0.530847 0.329338 61.186076 Choloepus_didactylus 0.537227 0.329232 63.175815 Vicugna_pacos 0.553177 0.338170 63.579561 Choloepus_hoffmanni 0.539357 0.329658 63.611076 Tursiops_truncatus 0.533317 0.325632 63.779051 Myrmecophaga_tridactyla 0.589117 0.357048 64.996583 # also calculated with SSREV model: mkdir /cluster/data/hg38/bed/cactus447/4dSSREV cd /cluster/data/hg38/bed/cactus447/4dSSREV cp -p ../4d/4d.all.mfa . cp -p ../4d/tree-commas.nh . time /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin/phyloFit \ --EM --precision MED --msa-format FASTA --subst-mod SSREV \ --tree tree-commas.nh 4d.all.mfa # real 1406m51.018s mv phyloFit.mod all.mod grep BACK all.mod # BACKGROUND: 0.254978 0.245022 0.245022 0.254978 ######################################################################### # phyloP for 447-way (DONE - 2017-11-06 - Hiram) # # split SS files into 1M chunks, this business needs smaller files # to complete ssh ku mkdir /hive/data/genomes/hg38/bed/cactus447/consPhyloP cd /hive/data/genomes/hg38/bed/cactus447/consPhyloP mkdir ss run.split cd run.split # the annotated maf's are in: ../iRows/result/chr*.maf printf '#!/bin/csh -ef set c = $1 set MAF = /hive/data/genomes/hg38/bed/cactus447/iRows/result/$c.maf set WINDOWS = /hive/data/genomes/hg38/bed/cactus447/consPhyloP/ss/$c set NL = `grep -m 10 -c -v "^#" $MAF` if ( -s $2 ) then exit 0 endif if ( -s $2.running ) then exit 0 endif date >> $2.running rm -fr $WINDOWS mkdir -p $WINDOWS pushd $WINDOWS > /dev/null if ( $NL > 0 ) then /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin/msa_split \ $MAF -i MAF -o SS -r $WINDOWS/$c -w 1000000,0 -I 1000 -B 5000 endif popd > /dev/null date >> $2 rm -f $2.running ' > doSplit.csh chmod +x doSplit.csh # do the easy ones first to see some immediate results ls -1SL -r ../../iRows/result | grep chr | sed -e "s/.maf//;" > maf.list # this needs a {check out line+ $(root1.done)} test for verification: printf '#LOOP ./doSplit.csh $(root1) {check out line+ $(root1).done} #ENDLOOP ' > template gensub2 maf.list single template jobList # this was tricky, some of the jobs couldn't complete even with # 128g, had to complete a number of these manually, and they # take a long time, more than 12 hours for individual chromosomes para -ram=128g create jobList para try ... check ... push ... etc... # Completed: 79 of 101 jobs # Crashed: 22 jobs # CPU time in finished jobs: 8058s 134.29m 2.24h 0.09d 0.000 y # IO & Wait Time: 380s 6.34m 0.11h 0.00d 0.000 y # Average job time: 107s 1.78m 0.03h 0.00d # Longest finished job: 3239s 53.98m 0.90h 0.04d # Submission to last job: 5350s 89.17m 1.49h 0.06d # run phyloP with score=LRT ssh ku mkdir /cluster/data/hg38/bed/cactus447/consPhyloP cd /cluster/data/hg38/bed/cactus447/consPhyloP mkdir run.phyloP cd run.phyloP # Adjust model file base composition background and rate matrix to be # representative of the chromosomes in play grep BACK ../../4dSSREV/all.mod # BACKGROUND: 0.254978 0.245022 0.245022 0.254978 # interesting, they are already paired up grep BACKGROUND ../../4dSSREV/all.mod | awk '{printf "%0.3f\n", $3 + $4}' # 0.490 /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin/modFreqs \ ../../4dSSREV/all.mod 0.490 > all.mod # verify, the BACKGROUND should now be paired up: grep BACK all.mod # BACKGROUND: 0.255000 0.245000 0.245000 0.255000 printf '#!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin set f = $1 set ssFile = $1:t set out = $2 set cName = $f:h set n = $f:r:e set grp = $cwd:t set cons = /hive/data/genomes/hg38/bed/cactus447/consPhyloP set tmp = $cons/tmp/$grp/$f /bin/rm -fr $tmp /bin/mkdir -p $tmp set ssSrc = "$cons/ss/$cName/$ssFile" set useGrp = "$grp.mod" /bin/ln -s $cons/run.phyloP/$grp.mod $tmp pushd $tmp > /dev/null echo source: $ssSrc.ss $PHASTBIN/phyloP --method LRT --mode CONACC --wig-scores --chrom $cName \ -i SS $useGrp $ssSrc.ss > $ssFile.wigFix popd > /dev/null /bin/mkdir -p $out:h sleep 4 /bin/touch $out:h /bin/mv $tmp/$ssFile.wigFix $out /bin/rm -fr $tmp /bin/rmdir --ignore-fail-on-non-empty $cons/tmp/$grp /bin/rmdir --ignore-fail-on-non-empty $cons/tmp ' > doPhyloP.csh chmod +x doPhyloP.csh # Create list of chunks find ../ss -type f | sed -e "s/.ss$//; s#../ss/##;" > ss.list # make sure the list looks good wc -l ss.list # 3031 ss.list # Create template file # file1 == $chr/$chunk/file name without .ss suffix printf '#LOOP ../run.phyloP/doPhyloP.csh $(path1) {check out line+ wigFix/$(dir1)/$(file1).wigFix} #ENDLOOP ' > template ###################### Running all species ####################### # setup run for all species mkdir /hive/data/genomes/hg38/bed/cactus447/consPhyloP/all cd /hive/data/genomes/hg38/bed/cactus447/consPhyloP/all mkdir wigFix gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList # beware overloading the cluster with these quick and high I/O jobs para -ram=64g create jobList para try ... check ... para -maxJob=16 push para time > run.time Completed: 3031 of 3031 jobs CPU time in finished jobs: 28585172s 476419.53m 7940.33h 330.85d 0.906 y IO & Wait Time: 132590s 2209.83m 36.83h 1.53d 0.004 y Average job time: 9475s 157.91m 2.63h 0.11d Longest finished job: 13733s 228.88m 3.81h 0.16d Submission to last job: 38568s 642.80m 10.71h 0.45d mkdir downloads time for D in `ls -d wigFix/chr* | sed -e 's#wigFix/##'` do echo "working: $D" 1>&2 find ./wigFix/${D} -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \ | sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \ | gzip -c > downloads/${D}.phyloP447way.wigFix.gz done # real 34m31.522s du -hsc downloads # 5.4G downloads # check integrity of data with wigToBigWig zcat downloads/*.wigFix.gz > all.wigFix time (wigToBigWig -verbose=2 all.wigFix /hive/data/genomes/hg38/chrom.sizes \ phyloP447way.bw) > bigWig.log 2>&1 egrep "real|VmPeak" bigWig.log # pid=236037: VmPeak: 33037308 kB # real 25m42.759s bigWigInfo phyloP447way.bw | sed -e 's/^/# /;' # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 7,824,240,203 # primaryIndexSize: 90,693,508 # zoomLevels: 10 # chromCount: 98 # basesCovered: 2,874,168,445 # mean: -0.022678 # min: -20.000000 # max: 8.123000 # std: 1.692685 # encode those files into wiggle data # time (zcat downloads/*.wigFix.gz \ wigEncode all.wigFix phyloP447way.wig phyloP447way.wib # Converted all.wigFix, upper limit 8.12, lower limit -20.00 # real 8m8.231s # -rw-rw-r-- 1 2874168445 Aug 30 14:08 phyloP447way.wib # -rw-rw-r-- 1 298978228 Aug 30 14:08 phyloP447way.wig du -hsc *.wi? # 2.7G phyloP447way.wib # 286M phyloP447way.wig # Load gbdb and database with wiggle. rm -f /gbdb/hg38/cactus447way/phyloP447way.wib ln -s `pwd`/phyloP447way.wib /gbdb/hg38/cactus447way/phyloP447way.wib time hgLoadWiggle -pathPrefix=/gbdb/hg38/cactus447way hg38 \ phyloP447way phyloP447way.wig # real 0m27.533s # use to set trackDb.ra entries for wiggle min and max # and verify table is loaded correctly wigTableStats.sh hg38 phyloP447way # db.table min max mean count sumData hg38.phyloP447way -20 8.123 -0.022678 2874168445 -6.51804e+07 # 1.69269 viewLimits=-8.4861:8.123 # stdDev viewLimits # that range is: 20+8.123= 28.123 for hBinSize=0.028123 # Create histogram to get an overview of all the data time hgWiggle -doHistogram \ -hBinSize=0.028123 -hBinCount=1000 -hMinVal=-20 -verbose=2 \ -db=hg38 phyloP447way > histogram.data 2>&1 # real 1m39.861s # x-axis range: grep -v chrom histogram.data | grep "^[0-9]" | ave -col=2 stdin \ | sed -e 's/^/# /;' # Q1 -10.086650 # median -4.012075 # Q3 2.062495 # average -4.017054 # min -20.000000 # max 8.123000 # count 864 # total -3470.734684 # standard deviation 7.028195 # find out the range for the 2:5 graph grep -v chrom histogram.data | grep "^[0-9]" | ave -col=5 stdin \ | sed -e 's/^/# /;' # Q1 0.000016 # median 0.000074 # Q3 0.000431 # average 0.001157 # min 0.000000 # max 0.025240 # count 864 # total 1.000001 # standard deviation 0.003283 # create plot of histogram: printf 'set terminal pngcairo size 1200,600 background "#000000" font "/usr/share/fonts/default/Type1/n022004l.pfb" set output "hg38.phyloP447.histo.png" set size 1.0, 1.0 set style line 1 lt 2 lc rgb "#ff88ff" lw 2 set style line 2 lt 2 lc rgb "#66ff66" lw 2 set style line 3 lt 2 lc rgb "#ffff00" lw 2 set style line 4 lt 2 lc rgb "#ffffff" lw 2 set border lc rgb "#ffff00" set key left box ls 3 set key tc variable set grid noxtics set grid ytics ls 4 set y2tics border nomirror out tc rgb "#ffffff" set ytics border nomirror out tc rgb "#ffffff" set title " Human hg38 Histogram phyloP447way track" tc rgb "#ffffff" set xlabel " hg38.phyloP447way score" tc rgb "#ffffff" set ylabel " Relative Frequency" tc rgb "#ff88ff" set y2label " Cumulative Relative Frequency (CRF)" tc rgb "#66ff66" set xrange [-4:2.5] set yrange [0:0.04] set y2range [0:1] plot "histogram.data" using 2:5 title " RelFreq" with impulses ls 1, \ "histogram.data" using 2:7 axes x1y2 title " CRF" with lines ls 2 ' | gnuplot # verify it looks sane display hg38.phyloP447.histo.png & # it now shows the spike at -1.000 for the artifical filling of # no data available. ###################### Running primate species ####################### # computing 4d-sites phyloFit for primates mkdir /cluster/data/hg38/bed/cactus447/4d/primates cd /cluster/data/hg38/bed/cactus447/4d/primates # the primates.list.txt was obtained from the hg38.447way.nh # tree by taking the first 243 in that tree, ending in: head -243 ../../hg38.447way.nh | tail -1 Galagoides_demidoff:0.00802126):0.0124128):0.0334321):0.0146219):0.105, head -243 ../../hg38.447way.nh | sed -e 's/^ \+//;' \ | tr -d '(' | sed -e 's/:.*//;' | sort > primates.list.txt # then only using the mfa sequences for those specific sequences ls ../mfa/chr*.mfa | while read C do B=`basename $C` faSomeRecords "${C}" primates.list.txt ${B} printf "%s\n" "${B}" done # resulting in: faSize chr*.mfa # 151201516 bases (18902 N's 151182614 real 151182614 upper 0 lower) # in 5832 sequences in 24 files # then running phyloFit: sed 's/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' primates.nh > tree-commas.nh # note: using the SSREV model here: /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin/phyloFit \ --EM --precision MED --msa-format FASTA --subst-mod SSREV \ --tree tree-commas.nh primates.mfa # rename the result mv phyloFit.mod primates.mod # adjust the frequencies cd /cluster/data/hg38/bed/cactus447/consPhyloP/run.phyloP grep BACK ../../4d/primates/primates.mod # BACKGROUND: 0.261935 0.238065 0.238065 0.261935 # interesting, they are already paired up grep BACK ../../4d/primates/primates.mod | awk '{printf "%0.3f\n", $3 + $4}' # 0.476 /cluster/bin/phast.build/cornellCVS/phast.2021-06-01/bin/modFreqs \ ../../4d/primates/primates.mod 0.476 > primates.mod # verify, the BACKGROUND should now be paired up: grep BACK all.mod # BACKGROUND: 0.262000 0.238000 0.238000 0.262000 # setup run for primate species mkdir /hive/data/genomes/hg38/bed/cactus447/consPhyloP/primates cd /hive/data/genomes/hg38/bed/cactus447/consPhyloP/primates mkdir wigFix gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList para -ram=6g create jobList para try ... check ... para -maxJob=16 push para time > run.time Completed: 2438 of 2438 jobs CPU time in finished jobs: 18916928s 315282.13m 5254.70h 218.95d 0.600 y IO & Wait Time: 68909s 1148.49m 19.14h 0.80d 0.002 y Average job time: 7787s 129.79m 2.16h 0.09d Longest finished job: 12050s 200.83m 3.35h 0.14d Submission to last job: 26707s 445.12m 7.42h 0.31d mkdir downloads time for D in `ls -d wigFix/chr* | sed -e 's#wigFix/##'` do echo "working: $D" 1>&2 find ./wigFix/${D} -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \ | sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \ | gzip -c > downloads/${D}.phyloP447wayPrimates.wigFix.gz done # real 33m18.006s du -hsc downloads # 4.9G downloads # check integrity of data with wigToBigWig zcat downloads/*.wigFix.gz > primates.wigFix time (wigToBigWig -verbose=2 primates.wigFix \ /hive/data/genomes/hg38/chrom.sizes \ phyloP447wayPrimates.bw) > bigWig.log 2>&1 egrep "real|VmPeak" bigWig.log # pid=236355: VmPeak: 33037308 kB # real 26m17.042s bigWigInfo phyloP447wayPrimates.bw | sed -e 's/^/# /;' # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 6,954,238,655 # primaryIndexSize: 90,693,508 # zoomLevels: 10 # chromCount: 98 # basesCovered: 2,874,168,445 # mean: -0.066955 # min: -20.000000 # max: 1.587000 # std: 1.479093 # encode those files into wiggle data # time (zcat downloads/*.wigFix.gz \ time wigEncode primates.wigFix phyloP447wayPrimates.wig phyloP447wayPrimates.wib # Converted primates.wigFix, upper limit 1.59, lower limit -20.00 # real 9m17.113s -XXX - running - Wed Aug 30 13:59:09 PDT 2023 # -rw-rw-r-- 1 2874168445 Aug 30 14:08 phyloP447wayPrimates.wib # -rw-rw-r-- 1 320006223 Aug 30 14:08 phyloP447wayPrimates.wig du -hsc *.wi? # 2.7G phyloP447way.wib # 306M phyloP447way.wig # Load gbdb and database with wiggle. rm -f /gbdb/hg38/cactus447way/phyloP447wayPrimates.wib ln -s `pwd`/phyloP447wayPrimates.wib /gbdb/hg38/cactus447way/phyloP447wayPrimates.wib time hgLoadWiggle -pathPrefix=/gbdb/hg38/cactus447way hg38 \ phyloP447wayPrimates phyloP447wayPrimates.wig # real 0m27.140s # use to set trackDb.ra entries for wiggle min and max # and verify table is loaded correctly wigTableStats.sh hg38 phyloP447wayPrimates # db.table min max mean count sumData hg38.phyloP447wayPrimates -20 1.587 -0.0669553 2874168445 -1.92441e+08 # 1.47909 viewLimits=-7.46242:1.587 # stdDev viewLimits # that range is: 20+1.587= 21.587 for hBinSize=0.021587 # Create histogram to get an overview of all the data time hgWiggle -doHistogram \ -hBinSize=0.021587 -hBinCount=1000 -hMinVal=-20 -verbose=2 \ -db=hg38 phyloP447wayPrimates > histogram.data 2>&1 # real 1m30.401s # x-axis range: grep -v chrom histogram.data | grep "^[0-9]" | ave -col=2 stdin \ | sed -e 's/^/# /;' # Q1 -11.645800 # median -7.231290 # Q3 -2.816750 # average -7.238782 # min -20.000000 # max 1.565410 # count 818 # total -5921.323540 # standard deviation 5.110470 # find out the range for the 2:5 graph grep -v chrom histogram.data | grep "^[0-9]" | ave -col=5 stdin \ | sed -e 's/^/# /;' # Q1 0.000007 # median 0.000032 # Q3 0.000340 # average 0.001237 # min 0.000001 # max 0.025761 # count 818 # total 1.012033 # standard deviation 0.003386 # create plot of histogram: printf 'set terminal pngcairo size 1200,600 background "#000000" font "/usr/share/fonts/default/Type1/n022004l.pfb" set output "hg38.phyloP447Primates.histo.png" set size 1.0, 1.0 set style line 1 lt 2 lc rgb "#ff88ff" lw 2 set style line 2 lt 2 lc rgb "#66ff66" lw 2 set style line 3 lt 2 lc rgb "#ffff00" lw 2 set style line 4 lt 2 lc rgb "#ffffff" lw 2 set border lc rgb "#ffff00" set key left box ls 3 set key tc variable set grid noxtics set grid ytics ls 4 set y2tics border nomirror out tc rgb "#ffffff" set ytics border nomirror out tc rgb "#ffffff" set title " Human hg38 Histogram phyloP447wayPrimates track" tc rgb "#ffffff" set xlabel " hg38.phyloP447wayPrimates score" tc rgb "#ffffff" set ylabel " Relative Frequency" tc rgb "#ff88ff" set y2label " Cumulative Relative Frequency (CRF)" tc rgb "#66ff66" set xrange [-3:2] set yrange [0:0.04] set y2range [0:1] plot "histogram.data" using 2:5 title " RelFreq" with impulses ls 1, \ "histogram.data" using 2:7 axes x1y2 title " CRF" with lines ls 2 ' | gnuplot # verify it looks sane display hg38.phyloP447Primates.histo.png & # it now shows the spike at -1.000 for the artifical filling of # no data available. #############################################################################