#For every run, we shall loop through all directories and run these commands /32 threads

#cd to specific thread dir per sample 
#For every dir touch comp_res.txt

#Sorting and creating indexes for bam files missing these
samtools sort A_sorted.keep.bam -o A_sorted.keep.bam
samtools index A_sorted.keep.bam A_sorted.keep.bai

samtools sort A_sorted.to.remap.bam -o A_sorted.to.remap.bam
samtools  index  A_sorted.to.remap.bam A_sorted.to.remap.bai 

#Collecting all vW tagged reads from STAR alignment #Note that not all reads have tags - only reads that overlap variants have tags and need remapping
samtools view A_sorted.bam | grep vW:i> STAR_vW_Tagged_Reads 
echo "STAR_vW_Tagged_Reads:" >> comp_res.txt #Writing out all stats to file
wc -l STAR_vW_Tagged_Reads >> comp_res.txt #None unique reads
echo " " >> comp_res.txt

#Extracting unique reads
samtools view A_sorted.bam | grep vW:i | sort -u -k1,1> STAR_vW_Tagged_Reads_Unique 
echo "STAR_vW_Tagged_Reads_Unique:" >> comp_res.txt
wc -l STAR_vW_Tagged_Reads_Unique  >> comp_res.txt
echo " " >> comp_res.txt

#Extracting only unique read-tag pairs
samtools view A_sorted.bam | grep vW:i | awk '{print $1 "\t" $(NF)}' | sort -u>  STAR_vW_Tagged_Reads_Unique_Subset

#Extracting other columns of interest
samtools view A_sorted.bam | grep vW:i | awk '{print $1 " " substr($0, index($0,$12)) }' | sort -u -k1,1>  STAR_vW_Tagged_Reads_Unique_Subset2 #Contains more attributes that we can leverage to further analyze the reads - here we are interested in the read id (column 1) and attribute columns, 12 to the last/NF 


#WASP reads to remap
samtools view A_sorted.to.remap.bam | awk 'BEGIN {while (getline < "STAR_vW_Tagged_Reads") w[$1]=$0; print length(w)} {if ($1 in w) {ww[$1]=0} else {print > "STAR_vW_Tagged_Reads_to_remap_novW"} } END {for (ii in w) {if (!(ii in ww)) print w[ii] > "STAR_vW_Tagged_Reads_notRemap"}}'
#Note that these are WASP reads that need remapping 

#Checking WASP reads that need remapping:
samtools view A_sorted.to.remap.bam | grep vW:i | awk '{print $1}' | sort -u> reads_to_remap_unique
echo "Unique reads to remap:" >> comp_res.txt
wc -l reads_to_remap_unique >> comp_res.txt #These are WASP reads that need remapping -  reads found to overlap variants in WASP but not in STAR
echo " " >> comp_res.txt

# STAR_vW_Tagged_Reads_to_remap_novW =  reads that were found to overlap variants by WASP but not STAR
# STAR_vW_Tagged_Reads_notRemap = vice versa
# STAR_vW_Tagged_Reads_to_remap_novW should not exist, i.e. all A.to.remap.bam reads have vW tags -  no system file ("STAR_vW_Tagged_Reads_to_remap_novW") generated either
echo "Number of reads overlapping variants by STAR but not WASP:" >> comp_res.txt
wc -l STAR_vW_Tagged_Reads_notRemap >> comp_res.txt # STAR_vW_Tagged_Reads_notRemap --> # number of reads overlapping variants by STAR but not WASP
echo " " >> comp_res.txt

#WASP reads to keep
samtools view A_sorted.keep.bam | grep vW:i | awk '{print $1}' | sort -u> reads_to_keep_unique
echo "Unique WASP reads to keep:" >> comp_res.txt
wc -l reads_to_keep_unique >> comp_res.txt
echo " " >> comp_res.txt

#Comparative matrix files: STAR_vW_Tagged_Reads_Unique_Subset(unique reads with tags), reads_to_remap_unique, reads_to_keep_unique
#Checking duplicates in all input files"
awk -F, 'a[$1]++{print $1}' STAR_vW_Tagged_Reads_Unique_Subset
awk -F, 'a[$1]++{print $1}' reads_to_remap_unique
awk -F, 'a[$1]++{print $1}' reads_to_keep_unique
#All files have no duplicate reads - Consider including branching logic for instances where we may have or need to consider duplications

awk 'FNR==NR{a[$1]=$1;next}{print $0,a[$1]?a[$1]:"NA"}' reads_to_remap_unique STAR_vW_Tagged_Reads_Unique_Subset2 > pass1_to_remap.txt # we can either tag reads that do not intersect as NA or 0: {print $0,a[$1]?a[$1]:0}
echo "Overlap between reads to remap and vW tagged reads" >> comp_res.txt
wc -l pass1_to_remap.txt >> comp_res.txt
echo " " >> comp_res.txt
#The output of pass1_to_remap.txt is  union join of vW_Tagged reads and reads that need remapping. The file is filled with NAs where a match is not found. e.g
#ERR1050081.10000004	vW:i:1 NA
#ERR1050081.10000031	vW:i:1 NA
#ERR1050081.10000048	vW:i:5 NA
#ERR1050081.10000085	vW:i:1 NA
#ERR1050081.10000105	vW:i:1 NA
#ERR1050081.10000238	vW:i:1 NA
#ERR1050081.10000311	vW:i:1 NA
#ERR1050081.10000460	vW:i:1 NA
#ERR1050081.10000491	vW:i:1 NA
#ERR1050081.10000525	vW:i:1 NA
#ERR1050081.10000576	vW:i:1 NA
#ERR1050081.1000058	vW:i:1 NA
#ERR1050081.10000595	vW:i:1 NA
#ERR1050081.10000613	vW:i:1 NA
#ERR1050081.10000634	vW:i:1 ERR1050081.10000634
#ERR1050081.10000667	vW:i:1 ERR1050081.10000667
#ERR1050081.10000676	vW:i:1 ERR1050081.10000676
#ERR1050081.10000678	vW:i:1 ERR1050081.10000678
#ERR1050081.10000736	vW:i:1 ERR1050081.10000736
#ERR1050081.1000078	vW:i:1 ERR1050081.1000078


awk 'FNR==NR{a[$1]=$1;next}{print $0,a[$1]?a[$1]:"NA"}' reads_to_keep_unique pass1_to_remap.txt > pass2_to_keep.txt
echo "Overlap between reads to keep and vW tagged reads:" >> comp_res.txt
wc -l pass2_to_keep.txt >> comp_res.txt
echo " " >> comp_res.txt
#The output of pass2_to_keep.txt is  union join of pass1_to_remap and reads that need keeping. Similarily colum 3 that represents reads to keep is filled with NAs where a match with vW_Tagged reads (column 1) is not found. e.g
#vW_Tagged reads	#Tag	#to remap 	#to keep
#ERR1050081.10000004	vW:i:1 NA ERR1050081.10000004
#ERR1050081.10000031	vW:i:1 NA ERR1050081.10000031
#ERR1050081.10000048	vW:i:5 NA NA
#ERR1050081.10000085	vW:i:1 NA ERR1050081.10000085
#ERR1050081.10000105	vW:i:1 NA ERR1050081.10000105
#ERR1050081.10000238	vW:i:1 NA ERR1050081.10000238
#ERR1050081.10000311	vW:i:1 NA NA
#ERR1050081.10000460	vW:i:1 NA ERR1050081.10000460
#ERR1050081.10000491	vW:i:1 NA ERR1050081.10000491
#ERR1050081.10000525	vW:i:1 NA ERR1050081.10000525
#ERR1050081.10000576	vW:i:1 NA NA
#ERR1050081.1000058	vW:i:1 NA ERR1050081.1000058
#ERR1050081.10000595	vW:i:1 NA ERR1050081.10000595
#ERR1050081.10000613	vW:i:1 NA NA
#ERR1050081.10000634	vW:i:1 ERR1050081.10000634 ERR1050081.10000634
#ERR1050081.10000667	vW:i:1 ERR1050081.10000667 ERR1050081.10000667
#ERR1050081.10000676	vW:i:1 ERR1050081.10000676 ERR1050081.10000676
#ERR1050081.10000678	vW:i:1 ERR1050081.10000678 NA
#ERR1050081.10000736	vW:i:1 ERR1050081.10000736 NA
#ERR1050081.1000078	vW:i:1 ERR1050081.1000078 NA