#sort your vcf file to match the headers in the reference genome
perl /projects/mkreitzman_prj/RNA-seq_simulation/GATK/vcfsorter.pl /projects/mkreitzman_prj/RNA-seq_simulation/GATK/Homo_sapiens_assembly19.dict /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.vcf > /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.vcf 2> /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/vcf_sort.err

# if you want to make a different reference for each haplotype, split the vcf into variants that are hetero- and homo- zygous
# make sure the PL field of the vcf != -1 (this causes GATK to crash)
perl /projects/mkreitzman_prj/RNA-seq_simulation/GATK/vcf_splitter.pl /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.vcf /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.homozygous.vcf /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.heterozygous.vcf

# use GATK to make a new reference for each haplotype using 1) the original vcf 2) the homozygous vcf only 
/gsc/software/linux-x86_64-centos5/java-1.7.0-u13/bin/java -jar /projects/analysis/analysis4/vcall/gatk/GenomeAnalysisTK.jar -R /home/pubseq/genomes/9606/hg19a/bwa_ind/genome/GRCh37-lite.fa -T FastaAlternateReferenceMaker -o /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1.fa --variant /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.vcf
/gsc/software/linux-x86_64-centos5/java-1.7.0-u13/bin/java -jar /projects/analysis/analysis4/vcall/gatk/GenomeAnalysisTK.jar -R /home/pubseq/genomes/9606/hg19a/bwa_ind/genome/GRCh37-lite.fa -T FastaAlternateReferenceMaker -o /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2.fa --variant /projects/mkreitzman_prj/RNA-seq_simulation/Illumina_body_map/HCT20170/ERS025093_5_lanes_dupsFlagged.sorted.homozygous.vcf

# munge the header lines generated in the GATK back into the correct fasta header lines (from the original reference)
gawk 'FNR==NR{/^>/ ? m=1 : m =0 ; if (m==1){b++; a[b]=$1;next} else {next}} { /^>/  ? m=1 : m=0 ; if (m==1) {c++; print a[c]} else {print $0}}' /home/pubseq/genomes/9606/hg19a/bwa_ind/genome/GRCh37-lite.fa /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1.fa >/projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1.with_headers.fa
rm /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1.fa
gawk 'FNR==NR{/^>/ ? m=1 : m =0 ; if (m==1){b++; a[b]=$1;next} else {next}} { /^>/  ? m=1 : m=0 ; if (m==1) {c++; print a[c]} else {print $0}}' /home/pubseq/genomes/9606/hg19a/bwa_ind/genome/GRCh37-lite.fa /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2.fa >/projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2.with_headers.fa
rm /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2.fa

# split the new reference into files (one chromosome per file) and add "chr" to each of the filenames
cd /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1_split_ref
awk '/^>/ {OUT=substr($1,2) ".fa"}; OUT {print >OUT}' /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap1.with_headers.fa
for i in *; do mv "$i" `echo "$i" | sed 's/^/chr/'`; done
cd /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2_split_ref
awk '/^>/ {OUT=substr($1,2) ".fa"}; OUT {print >OUT}' /projects/mkreitzman_prj/RNA-seq_simulation/GATK/alt_ref_hap2.with_headers.fa
for i in *; do mv "$i" `echo "$i" | sed 's/^/chr/'`; done

# run flux simulator on first haplotype
/home/mkreitzman/Programs/flux-simulator-1.2.1/bin/flux-simulator -t simulator -x -l -s -p /projects/mkreitzman_prj/RNA-seq_simulation/Flux_simulator/hap1.par

# run flux simulator on second haplotype using only -s and other reference genome
/home/mkreitzman/Programs/flux-simulator-1.2.1/bin/flux-simulator -s -p /projects/mkreitzman_prj/RNA-seq_simulation/Flux_simulator/hap2.par

#rename the second haplotype  to avoid having 2 fragments with the smae read name
cat hap2.fastq | awk '{if ($0~/^@/) {gsub(/^@/,"@:hap2:"); print} else {print $0}}' >  hap2.renamed.fastq

#merge the fastqs for the two haplotypes
cd /projects/mkreitzman_prj/RNA-seq_simulation/Flux_simulator
cat hap1.fastq hap2.renamed.fastq > all.fastq

#separate merged fastq into ist and 2nd end files. 
perl /projects/mkreitzman_prj/RNA-seq_simulation/split_interleaved_fastq.pl
