GitLab

• Rather than putting the following in the "analysis" snakefile, can we code this into the Rscript when making the figures?

# using the mappability_index files. example: `hybrid_mmi_merged.txt` file in the `results/mapping` folder, which will be created using the snakefile I'm making.

for file in *.txt; do echo $file; egrep 'total | mapped' $file | head -n 2; done | paste - - - | \
   sed $'1 i\\\nsample\ttotal\tmapped\tpercent_mapped' | sed 's/ + 0 in total (QC-passed reads + QC-failed reads)//g' | \
     sed s'/ + 0 mapped//g' | sed 's/ : N\/A)//g' | sed s'/(//g' | sed s'/.txt//g' | \
       awk '{print$0=$1"\t"$2"\t"$3"\t"$4}' > mappability_index.tsv

And the same for the counts where contigs smaller than 1Kb - 5kb are removed. examples below

for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 1000 {sum += $3} END {print sum}' "$file"  # sum of number of reads (3rd column) if 2nd column is greater/equal to 1000
  awk '{sum+=$3+$4} END {print sum}' "$file"            # sum of 3rd and 4th columns to get total reads
done | paste - - - | awk '$4=100*$2/$3' | \             # getting the percentage from columns 2 and 3 and writing to column-4
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_1000bp.txt

# collating the sums for contigs larger than 1500 bp
for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 1500 {sum += $3} END {print sum}' "$file"
  awk '{sum+=$3+$4} END {print sum}' "$file"
done | paste - - - | awk '$4=100*$2/$3' | \
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_1500bp.txt

@valentina.galata @cedric.laczny what do you think? or should I try to code this into the snakefile? The "headache" is that the names of the BAM files are super long and different. Could create symlinks like I did in the past, but want to avoid creating more than we need to.