39b70f7f9425b87eb06dc521f2579abf86d877b6 hiram Wed Apr 12 09:13:06 2023 -0700 calculating coverage on hg38 of HPRC lift.over chains refs #30508 diff --git src/hg/makeDb/doc/hg38/hprcCoverage.txt src/hg/makeDb/doc/hg38/hprcCoverage.txt new file mode 100644 index 0000000..dd611d0 --- /dev/null +++ src/hg/makeDb/doc/hg38/hprcCoverage.txt @@ -0,0 +1,196 @@ +## calculate coverage measurement from liftOver chain files + +## collect files together in one location + +mkdir /hive/data/genomes/hg38/bed/hprc/liftOver/chains +cd /hive/data/genomes/hg38/bed/hprc/liftOver/chains + +cut -f1 ~/kent/src/hg/makeDb/doc/hprcAsmHub/hprc.orderList.tsv \ + | cut -d'_' -f1-2 > accessionId.txt + +# need to do the 'realpath' to the /hive/data/ location to be available +# for cluster runs. /gbdb/hg38 isn't available for cluster reference + +for id in `cat accessionId.txt` +do + p0="/gbdb/hg38/liftOver/hg38To${id}.over.chain.gz" + ln -s `realpath "${p0}"` ./ +done + +# and one extra for the T2T + +ln -s /hive/data/genomes/hg38/bed/lastzGCA_009914755.4.2022-03-28/axtChain/hg38.GCA_009914755.4.over.chain.gz \ + hg38ToGCA_009914755.4.over.chain.gz + +# should be 97: + +ls *.over.chain.gz | wc -l + 97 + +## chrom.sizes were previously obtained by Max: + +mkdir /hive/data/genomes/hg38/bed/hprc/sizes +cd /hive/data/genomes/hg38/bed/hprc/sizes +for i in `cat ../accs.txt`; do wget `curl -s "https://genome-test.gi.ucsc.edu/list/files?genome=$i&format=text" | grep sizes`; done + +## and the T2T +cp -p /hive/data/genomes/asmHubs/genbankBuild/GCA/009/914/755/GCA_009914755.4_T2T-CHM13v2.0/GCA_009914755.4_T2T-CHM13v2.0.chrom.sizes \ + GCA_009914755.4.chrom.sizes.txt + +## and the 2bit files + +mkdir /hive/data/genomes/hg38/bed/hprc/2bits +cd /hive/data/genomes/hg38/bed/hprc/2bits + +cut -f1 ~/kent/src/hg/makeDb/doc/hprcAsmHub/hprc.orderList.tsv \ + | cut -d'_' -f1-2 > accessionId.txt + +for id in `cat accessionId.txt` +do + export id + for P in `curl -s -L "https://genome-test.gi.ucsc.edu/list/files?genome=$id;format=text" \ + | grep 2bit | grep -v bpt | sed -e 's#https://hgdownload.soe.ucsc.edu/hubs#/hive/data/genomes/asmHubs##;'` + do + rp=`realpath $P` + echo ln -s \"${rp}\" ./$id.2bit + ln -s "${rp}" ./$id.2bit + done +done + +## and the T2T 2bit file: +cp -p /hive/data/genomes/asmHubs/genbankBuild/GCA/009/914/755/GCA_009914755.4_T2T-CHM13v2.0/GCA_009914755.4_T2T-CHM13v2.0.2bit \ + GCA_009914755.4.2bit + +## convert the chain files to bed: + +mkdir /hive/data/genomes/hg38/bed/hprc/chainToBed +cd /hive/data/genomes/hg38/bed/hprc/chainToBed + +## make a cluster jobList +ls ../liftOver/chains | grep chain.gz | while read C +do + s="${C#hg38To}" + c="${s%.over.chain.gz}" + printf "./runOne %s {check out exists+ %s.bed}\n" "${c}" "${c}" +done > jobList + +printf '#!/bin/bash + +set -beEu -o pipefail + +export S="${1}" +export overChain="/hive/data/genomes/hg38/bed/hprc/liftOver/chains/hg38To${S}.over.chain.gz" +export twoBit="/hive/data/genomes/hg38/bed/hprc/2bits/${S}.2bit" +export sizes="/hive/data/genomes/hg38/bed/hprc/sizes/${S}.chrom.sizes.txt" + +if [ "${overChain}" -nt "${S}.bed" ]; then + /cluster/bin/x86_64/chainToPsl $overChain /hive/data/genomes/hg38/chrom.sizes \ + $sizes /hive/data/genomes/hg38/hg38.2bit $twoBit stdout \ + | /cluster/bin/x86_64/pslToBed stdin stdout \ + | /cluster/bin/x86_64/bedToExons stdin ${S}.bed + touch -r $overChain ${S}.bed +fi +' > runOne +chmod +x runOne + +## running on ku +ssh ku +cd /hive/data/genomes/hg38/bed/hprc/chainToBed +para create jobList +para push +para time +Completed: 96 of 96 jobs +CPU time in finished jobs: 54194s 903.24m 15.05h 0.63d 0.002 y +IO & Wait Time: 231s 3.84m 0.06h 0.00d 0.000 y +Average job time: 567s 9.45m 0.16h 0.01d +Longest finished job: 1378s 22.97m 0.38h 0.02d +Submission to last job: 2286s 38.10m 0.64h 0.03d + +## and the T2T lift.over: +ssh hgwdev +cd /hive/data/genomes/hg38/bed/hprc/chainToBed + +./runOne GCA_009914755.4 + +## put all the bed files into one file +mkdir ../coverage +$HOME/bin/x86_64/gnusort -k1,1 -k2,2n -S100G --parallel=32 *.bed \ + > ../coverage/allOvers.bed + +## calculate the coverage numers +cd /hive/data/genomes/hg38/bed/hprc/coverage +bedItemOverlapCount hg38 allOvers.bed > allOvers.bedGraph +touch -r allOvers.sorted.bed allOvers.bedGraph + +## account for those areas not covered at all + +bedInvert.pl /hive/data/genomes/hg38/chrom.sizes allOvers.bedGraph \ + > zeroCounts.bed +touch -r allOvers.bedGraph zeroCounts.bed + +## add in the zero counts + +(awk '{printf "%s\t%d\t%d\t0\n", $1,$2,$3}' zeroCounts.bed; + cat allOvers.bedGraph) | $HOME/bin/x86_64/gnusort -k1,1 -k2,2n \ + -S100G --parallel=32 > fullGenome.bedGraph + +## final answer: + +bedGraphToBigWig fullGenome.bedGraph /hive/data/genomes/hg38/chrom.sizes \ + hprcCoverage.bw + +## count up percent of genome in each coverage number + +printf '#!/usr/bin/env perl + +use strict; +use warnings; + +my $genomeSize = `ave -col=2 /hive/data/genomes/hg38/chrom.sizes | grep -w total`; +chomp $genomeSize; +$genomeSize =~ s/.000000//; +$genomeSize =~ s/total +//; +printf "# hg38 total genome size: %s\n", $genomeSize; + +my $countedBases = 0; +my @genomeProfile; # index N is the coverage depth 0 to 97 + # value is total bases at that coverage +open (FH, "< fullGenome.bedGraph") or die "can not read fullGenome.bedGraph"; +while (my $line = <FH>) { + chomp $line; + my ($chr, $start, $end, $cover) = split('\s+', $line); + $countedBases += $end - $start; + $genomeProfile[$cover] += $end - $start; +} +close (FH); + +printf STDERR "# counted bases: %d\n", $countedBases; + +for (my $i = 0; $i < scalar(@genomeProfile); ++$i) { + printf "%d\t%d\t%.4f\n", $i, $genomeProfile[$i], 100.0 * $genomeProfile[$i] / $genomeSize; +} +' > percentGenome.pl + +chmod +x percentGenome.pl + +./percentGenome.pl > profile.txt + +## created an ordered by size bed file: + +awk '{printf "%s\t%d\t%d\t%d\t%d\n", $1,$2,$3,$4,$3-$2}' fullGenome.bedGraph \ + | sort -k5,5nr > bySize.bed5 + +## profile sizes of blocks for under 20 coverage + +awk '$4 < 20' bySize.bed5 | ave -col=5 stdin + +Q1 1.000000 +median 1.000000 +Q3 3.000000 +average 2358.696411 +min 1.000000 +max 30180026.000000 +count 97609 +total 230229998.000000 +standard deviation 172056.394886 +