genomewalker
10/7/2019 - 2:46 PM

Get PMD score

# Get list of contigs
bedtools genomecov -ibam ar3_10.2_gtdb_r89.bam -bg  \
| awk '$4 > 9{split($0,a,"\\."); print a[1]}' \
| sort -u > ar3_10.2_gtdb_r89.acc

# Get 
head -n10 ar3_10.2_gtdb_r89.acc \
| while read LINE; do 
  samtools view -b -L <(grep -w -F ${LINE} /vol/scratch/DBs/bowtie2/gtdb-r89_54k_genomes.fasta.fai \
  | awk '{printf("%s\t0\t%s\n",$1,$2);}') ar3_10.2_gtdb_r89.bam \
  | samtools view -h \
  | python pmdtools.0.60.py --writesamfield --threshold -100 --header \
  | samtools view -S \
  | awk '{ print $1"\t"$NF }' \
  | sed -e 's/DS:Z://g'; 
done > test_damage.txt