ml use $WORK/public/apps/modulefiles
ml repliscan
module load samtools bedtools/2.25.0 R
set -xe

REF=hg19_masked.fasta
if [ ! -e $REF ]
then
	# Download tar of chromosomes
	wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFaMasked.tar.gz
	for chrom in `tar -tnzf chromFaMasked.tar.gz | grep -v "_"`
	do
		tar -Oxzf chromFaMasked.tar.gz $chrom >> $REF
	done
	# Index fasta
	samtools faidx $REF
fi

# Download Hansen data
while read url
do
	F=$(basename $url)
	[ -e $F ] || curl -O -L $url
done << EOF
https://www.encodeproject.org/metadata/type=Experiment&assay_term_name=Repli-seq&replicates.library.biosample.donor.organism.scientific_name=Homo+sapiens&biosample_term_name=foreskin+fibroblast/metadata.tsv
https://www.encodeproject.org/files/ENCFF001GIH/@@download/ENCFF001GIH.bam
https://www.encodeproject.org/files/ENCFF001GIJ/@@download/ENCFF001GIJ.bam
https://www.encodeproject.org/files/ENCFF001GIT/@@download/ENCFF001GIT.bam
https://www.encodeproject.org/files/ENCFF001GIW/@@download/ENCFF001GIW.bam
https://www.encodeproject.org/files/ENCFF001GJF/@@download/ENCFF001GJF.bam
https://www.encodeproject.org/files/ENCFF001GJG/@@download/ENCFF001GJG.bam
https://www.encodeproject.org/files/ENCFF001GJA/@@download/ENCFF001GJA.bam
https://www.encodeproject.org/files/ENCFF001GJB/@@download/ENCFF001GJB.bam
https://www.encodeproject.org/files/ENCFF001GIA/@@download/ENCFF001GIA.bam
https://www.encodeproject.org/files/ENCFF001GIB/@@download/ENCFF001GIB.bam
https://www.encodeproject.org/files/ENCFF001GIN/@@download/ENCFF001GIN.bam
https://www.encodeproject.org/files/ENCFF001GIO/@@download/ENCFF001GIO.bam
EOF

# Rename based on metadata
grep "\.bam" metadata.tsv | cut -f 1,15,29 | awk '{print "for f in "$1"*; do ln -fs $f ${f/"$1"/"$2"_"$3"}; done"}' | xargs -n 1 -I {} bash -c {}

# Make 150bp windows of reference
export LC_ALL=C
BED=${REF}.150.bed
[ -e $BED ] || ( bedtools makewindows -g ${REF}.fai -w 150 | sort -S 20G -k1,1 -k2,2n > $BED )

# Combine replicates and remove any reads where the coverage in a 150bp window
# is greater than 4
for prefix in `ls [GS]*bam | sed -e "s/_[12].bam//" | sort -u`
do
	BG=${prefix}.bedgraph
	if [ ! -e $BG ]
	then
		for bam in ${prefix}*bam
		do
			bedtools bamtobed -i $bam
		done | sort -S 20G -k1,1 -k2,2n | bedtools intersect -a $BED -b - -c -sorted | awk '{if ($4 < 5) print $0;}' > tmp.bed
		mv tmp.bed $BG
	fi
done

#[ -e $SLIDE ] || bedtools makewindows -g ${REF}.fai -w 50000 -s 1000 | sort -S 1G -k1,1 -k2,2n > $SLIDE
# Make 1kb window reference
RBED=${REF}.1000.bed
[ -e $RBED ] || bedtools makewindows -g ${REF}.fai -w 1000 | sort -S 1G -k1,1 -k2,2n > $RBED

# Make bedgraphs normalized by percentage of total replication along the sliding windows
for bed in [GS]*[^e0]\.bedgraph
do
	OUT=${bed/.bed/_1000.bed}
	if [ ! -e $OUT ]
	then
		TOTAL=$(awk '{sum+=$4;} END {print sum;}' $bed)
		bedtools map -a $RBED -b $bed -o sum -c 4  | awk -v sum=$TOTAL '{
		printf("%s\t%i\t%i\t%f\n", $1, $2, $3, $4*100/sum);
	}' > tmp.bedgraph
	bioitools smooth -m average -b 50 tmp.bedgraph > $OUT
	rm tmp.bedgraph
	fi
done

###################################
# Run repliscan
###################################

INPUT=../input.txt
cat << EOF > input.txt
TS	G1b_1.bam	G1b_2.bam	G2_1.bam	G2_2.bam	S1_1.bam	S1_2.bam	S2_1.bam	S2_2.bam	S3_1.bam	S3_2.bam	S4_1.bam	S4_2.bam
G1b	G1b_1.bam	G1b_2.bam
S1	S1_1.bam	S1_2.bam
S2	S2_1.bam	S2_2.bam
S3	S3_1.bam	S3_2.bam
S4	S4_1.bam	S4_2.bam
G2	G2_1.bam	G2_2.bam
EOF

for bin in 50
do
	BED=${REF}.${bin}000.bed
	[ -e $BED ] || ( bedtools makewindows -g ${REF}.fai -w ${bin}000 | sort -S 5G -k1,1 -k2,2n > tmp.bed && mv tmp.bed $BED )
	fold=repliscan_${bin}kb
	[ -e $fold ] && rm -rf $fold
	mkdir $fold && ln -s $PWD/[GS]*bam $fold/
	cd $fold && repliscan.py -R norm -a sum -r ../$REF -w ${bin}000 --plot $INPUT &> log.out || echo $fold 'failed'
done
