#!/bin/tcsh

#Command pipeline for tallying polymorphic positions.

#These commands are assumed to be executed in a directory which contains subdirectories for each sample.
#In this example the path to this directory is ~/analysis
#Each subdirectory should have a name unique to each sample, which will be represented as $file
#Within each subdirectory are reads to be processed for that sample (here "$file.miss.fq").
#The commands fastq_to_fasta, fastq_quality_filter, and fastq_masker are part of the FASTX-Toolkit suite.
#fastq_collapse.py and polymorphic_read_counter_bwaPileup.pl are available at https://github.com/listonlab
#To personalize this analysis you will need to make changes on the following lines:
#Line 16: Change this range to match the range of the sequence you're interested in (in this example the reference sequence is 5839 bp).
#Line 17: Change the command within the parentheses to match all of your subdirectories (in this case all of the subdirectories began with 'Asclepias').
#Line 35: Change 2013 to the year the analysis is performed, this should match any directory that is created by alignreads.py
#Line 62: Change the inequalities (here 1 and 5839) to match the range of the sequence you're interested in

cd ~/analysis
seq 1 5839 | sort > base_proportions_concatenated.sort
foreach file ( `ls -1d --indicator-style=none Asclepias* ` )
cd ~/analysis/$file

###### Removes duplicate reads
fastq_collapse.py -i $file.miss.fq -o $file.miss.collapse.fq
fastq_to_fasta -n -v -i $file.miss.collapse.fq -o $file.miss.collapse.fa

###### Alignreads
python alignreads.py -w 5 -o linear -f alnrds.$file $file.miss.collapse.fa ../trim_syriaca5_cistron_ref.fsa

###### Quality filtering and masking
mkdir bwa
fastq_quality_filter -q 20 -p 01 -i $file.miss.collapse.fq -o $file.readqual_out.fq -v
fastq_masker -q 20 -i $file.readqual_out.fq -o $file.basequal_out.fq -v
mv $file.basequal_out.fq bwa/

###### Grabs the alignment between the group reference and the individual consensus
#alignreads.py created a directory with a name that contains the date of the analysis, in this example the commands were run in 2013, so all of these directories end in 2013
cd *2013
awk '/>Consensus/,/>Masked Consensus/' *.delta.qual.fsa | sed s/\>Masked\ Consensus// | sed s/\>Consensus/\>$file/ | sed s/-//g > $file.cistron_ref.fa
awk 'BEGIN {x=0} "'$5'"~/^[A-Z a-z]"'$'"/ {print "'$1'",++x,"'$4'","'$5'"}' *.delta.qual > $file.position2ref
mv $file.cistron_ref.fa ~/analysis/$file/bwa/
mv $file.position2ref ~/analysis/$file/bwa/

###### Read mapping with BWA and SAMtools
cd ~/analysis/$file/bwa
bwa index $file.cistron_ref.fa
bwa aln -f $file.sai $file.cistron_ref.fa $file.basequal_out.fq
bwa samse -f $file.sam $file.cistron_ref.fa $file.sai $file.basequal_out.fq
samtools view -bS $file.sam > $file.bam
samtools sort $file.bam $file.sort
samtools index $file.sort.bam
samtools mpileup -f $file.cistron_ref.fa  $file.sort.bam > $file.pileup.txt

###### Counting and parsing mapped reads with polymorphic_read_counter_bwaPileup.pl
polymorphic_read_counter_bwaPileup.pl $file.pileup.txt > $file.props
mkdir proportions
cd proportions
cut -f2,4 ../$file.props | sort -k1,1 > $file.props.sort
sort -k2,2 ../$file.position2ref > $file.position2ref.sort
join -1 2 -2 1 -o 1.2,1.3,1.4,2.2,1.1 --check-order $file.position2ref.sort $file.props.sort > $file.merged
sort -k5,5 $file.merged > $file.merged.sort
uniq -f4 $file.merged.sort > $file.merged.uniq
awk '{print $5"\t"$4}' $file.merged.uniq > $file.tojoin
#In this example we are interested in reference positions 1-5839. Change the ranges in the next command to match your range.
awk '{if (($1 >= 1) && ($1 <= 5839)) print $0}' $file.tojoin > $file.cistron.sort
sort -k1,1 -n $file.cistron.sort > $file.cistron

###### Adding the proportions at each position to a table containing information for each sample
cd ~/analysis
join -j1 --check-order -a1 base_proportions_concatenated.sort ./$file/bwa/proportions/$file.cistron.sort > props.temp
rm -f base_proportions_concatenated.sort
mv props.temp base_proportions_concatenated.sort

###### Closes the foreach loop
end
