diff --git a/2-generate_gaps/extract_gaps.py b/2-generate_gaps/extract_gaps.py index 416d547..06fff5a 100644 --- a/2-generate_gaps/extract_gaps.py +++ b/2-generate_gaps/extract_gaps.py @@ -271,7 +271,7 @@ def interval_merge(intervals): # splitting up the reference genome into the gaps ######################################################## -fasta_fn = open(sys.argv[2]) +fasta_fn = open(sys.argv[2], "r") output_fn = open(sys.argv[4], "w") gap_offset = 0 diff --git a/4-extract_reads/extract_reads.sh b/4-extract_reads/extract_reads.sh index 97dab70..69f49e2 100644 --- a/4-extract_reads/extract_reads.sh +++ b/4-extract_reads/extract_reads.sh @@ -4,5 +4,15 @@ READ_BAM_FILE=$1 BED_FILE=$2 READSIZE=$3 -awk -v bam_file="$READ_BAM_FILE" -v read_size="${READSIZE}M" '{cmd="samtools view -h "bam_file" " $1":"$2"-"$3" | convert2bed --input=sam -"; while(( cmd | getline result ) > 0 ){split(result, array, "\t"); if(array[2] >= ($2-1) && array[3] <= ($3-1) && (array[5] <= 10 || array[8] != read_size) ){if(and(array[7], 64)){print array[4]".1" "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}else if(and(array[7], 128)){print array[4]".2" "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}else{print array[4] "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}}} close(cmd);}' ${BED_FILE} | sort -uk1,1 +awk -v bam_file="$READ_BAM_FILE" -v read_size="${READSIZE}M" \ +'{cmd="samtools view -h "bam_file" " $1":"$2"-"$3" | convert2bed --input=sam -"; \ +while(( cmd | getline result ) > 0 ){\ + split(result, array, "\t"); \ + if(array[2] >= ($2-1) && array[3] <= ($3-1) && (array[5] <= 10 || array[8] != read_size) ){\ + if(and(array[7], 64)){print array[4]".1" "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}\ + else if(and(array[7], 128)){print array[4]".2" "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}\ + else{print array[4] "\t" array[2] "\t" array[3] "\t" array[1] "\t" array[5] "\t" array[8];}\ + }\ +} \ +close(cmd);}' ${BED_FILE} | sort -uk1,1 diff --git a/4-extract_reads/parallel_extract_reads.sh b/4-extract_reads/parallel_extract_reads.sh new file mode 100644 index 0000000..16f0759 --- /dev/null +++ b/4-extract_reads/parallel_extract_reads.sh @@ -0,0 +1,22 @@ +#!/bin/bash +# Need GNU parallel + +module load parallel +source $(which env_parallel.bash) + +#variables to change are: +#{BED_FILES_FOLDER} +#{READ_BAM_FILE} +#{EXTRACT_READS_SCRIPT} +#{NUM_FILES} +#{THREAD} + +FILES=({BED_FILES_FOLDER}/*.bed) +BED_FILE=${FILES[${SLURM_ARRAY_TASK_ID}]} +OUT_READS=$(echo ${BED_FILE}.reads) + +# Check the parallel command first +env_parallel -j3 -k --dryrun /usr/bin/time -v -p -o "${OUT_READS}.time" bash {EXTRACT_READS_SCRIPT} {READ_BAM_FILE} {1} {READSIZE} '>' {2} ::: $(echo ${FILES[@]}) ::: $(echo ${FILES[@]} | awk '{for(i=1;i<=NF;i++) printf $i".reads ";}') +env_parallel -j3 -k --progress /usr/bin/time -v -p -o "${OUT_READS}.time" bash {EXTRACT_READS_SCRIPT} {READ_BAM_FILE} {1} {READSIZE} '>' {2} ::: $(echo ${FILES[@]}) ::: $(echo ${FILES[@]} | awk '{for(i=1;i<=NF;i++) printf $i".reads ";}') + +module unload parallel diff --git a/4-extract_reads/slurm_extract_reads.sh b/4-extract_reads/slurm_extract_reads.sh index 62dc93e..bf5d4ff 100644 --- a/4-extract_reads/slurm_extract_reads.sh +++ b/4-extract_reads/slurm_extract_reads.sh @@ -18,6 +18,11 @@ FILES=({BED_FILES_FOLDER}/*.bed) BED_FILE=${FILES[${SLURM_ARRAY_TASK_ID}]} OUT_READS=$(echo ${BED_FILE}.reads) -srun --output=${OUT_READS} /usr/bin/time -v -p -o "${OUT_READS}.time" sh {EXTRACT_READS_SCRIPT} {READ_BAM_FILE} ${BED_FILE} {READSIZE} +srun --output=${OUT_READS} \ +/usr/bin/time -v -p -o "${OUT_READS}.time" \ +bash {EXTRACT_READS_SCRIPT} \ +{READ_BAM_FILE} \ +${BED_FILE} \ +{READSIZE} sleep 1 diff --git a/run_pipeline.sh b/run_pipeline.sh index a93c6e7..8d9b650 100644 --- a/run_pipeline.sh +++ b/run_pipeline.sh @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash SRCFOLDER=$1 #ref folders containing ref file: @@ -11,29 +11,298 @@ READ_BAM=$6 FIRST_PAIR=$7 SECOND_PAIR=$8 OUTPUT=$9 -THREAD=$10 - -mkdir -p "${OUTPUT}/" -#sh "${SRCFOLDER}/1-generate_chain/chain_install.sh" "${OUTPUT}/" -#sh "${SRCFOLDER}/1-generate_chain/sample_run.sh" "${OLDREF}" "${NEWREF}" "${OUTPUT}/" "${SRCFOLDER}/1-generate_chain/chain_generate.sh" ${SEQ_FILE_EXT} "${OUTPUT}/" -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/${chr}_extracted_gaps.fa.time" python "${SRCFOLDER}/2-generate_gaps/extract_gaps.py" $i "${NEWREF}/$chr.${SEQ_FILE_EXT}" ${READSIZE} "${OUTPUT}/${chr}_extracted_gaps.fa"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/${chr}_kmer_gaps.fasta.time" python "${SRCFOLDER}/2-generate_gaps/gaps_to_fasta.py" "${OUTPUT}/${chr}_extracted_gaps.fa" ${READSIZE} "${OUTPUT}/${chr}_kmer_gaps.fasta" 10; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; sh "${SRCFOLDER}/3-align_gaps/align_gaps.sh" "${OLDREF}/${chr}.${SEQ_FILE_EXT}" "${OUTPUT}/${chr}_kmer_gaps.fasta" ${THREAD} "${OUTPUT}/${chr}_aligned_gaps"; done -mkdir "${OUTPUT}/bedfiles"; -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; mkdir -p "${OUTPUT}/bedfiles/${chr}/"; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_regions.time" sh "${SRCFOLDER}/4-extract_reads/extract_regions.sh" "${OUTPUT}/${chr}_aligned_gaps.bam" "${OUTPUT}/bedfiles/${chr}/"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; num_of_files=`ls -l ${OUTPUT}/bedfiles/${chr}/*.bed | wc -l`; sed_bed_files=$(echo ${OUTPUT}\/bedfiles\/${chr}/ | sed 's/\//\\\//g'); sed_read_scripts=$(echo ${SRCFOLDER}\/4-extract_reads\/extract_reads.sh | sed 's/\//\\\//g'); sed_read_bam=$(echo $READ_BAM | sed 's/\//\\\//g'); cat "${SRCFOLDER}/4-extract_reads/slurm_extract_reads.sh" | sed "s/{NUM_FILES}/$num_of_files/g" | sed "s/{THREAD}/$THREAD/g" | sed "s/{BED_FILES_FOLDER}/$sed_bed_files/g" | sed "s/{EXTRACT_READS_SCRIPT}/$sed_read_scripts/g" | sed "s/{READ_BAM_FILE}/$sed_read_bam/g" | sed "s/{READSIZE}/$READSIZE/g" > "${OUTPUT}/${chr}_samplejob.sh"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; sbatch --wait "${OUTPUT}/${chr}_samplejob.sh"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed.time" cat "${OUTPUT}/bedfiles/${chr}/"*.bed | sortBed -i - | mergeBed -i - > "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; mkdir -p "${OUTPUT}/bedfiles/${chr}/retired_bed/"; mkdir -p "${OUTPUT}/bedfiles/${chr}/constant_bed/"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed.time" python "${SRCFOLDER}/4-extract_reads/get_retired_regions.py" ${READSIZE} "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed" $i "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed.time" sh "${SRCFOLDER}/4-extract_reads/extract_reads_noprune.sh" ${READ_BAM} "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed" > "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; mkdir -p "${OUTPUT}/bedfiles/${chr}/reads/"; cat "${OUTPUT}/bedfiles/${chr}/"*.reads "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed" > "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed"; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_sequences_full.time" sh "${SRCFOLDER}/4-extract_reads/extract_sequence.sh" ${FIRST_PAIR} ${SECOND_PAIR} "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed" "${OUTPUT}/bedfiles/${chr}/reads/"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.time" bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/reads_1.fastq" "${OUTPUT}/bedfiles/${chr}/reads/reads_2.fastq" | samtools view -h -F4 | samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.bam"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.time" bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" | samtools view -h -F4 | samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.bam"; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; sh "${SRCFOLDER}/0-align_reads.sh" "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/reads" "${OUTPUT}/bedfiles/${chr}/reads/paired" $THREAD; done -for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; sh "${SRCFOLDER}/0-align_singletons.sh" "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" "${OUTPUT}/bedfiles/${chr}/reads/singletons" $THREAD; done +THREAD=${10} + + +function generate_chain(){ + bash "${SRCFOLDER}/1-generate_chain/chain_install.sh" "${OUTPUT}" + bash "${SRCFOLDER}/1-generate_chain/sample_run.sh" \ + "${OLDREF}" \ + "${NEWREF}" \ + "${OUTPUT}/" \ + "${SRCFOLDER}/1-generate_chain/chain_generate.sh" \ + ${SEQ_FILE_EXT} \ + "${OUTPUT}/" +} + +function generate_gaps(){ + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//`; \ + /usr/bin/time -v -p -o "${OUTPUT}/${chr}_extracted_gaps.fa.time" \ + python3 "${SRCFOLDER}/2-generate_gaps/extract_gaps.py" \ + $i \ + "${NEWREF}/$chr.${SEQ_FILE_EXT}" \ + ${READSIZE} \ + "${OUTPUT}/${chr}_extracted_gaps.fa" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//`; \ + /usr/bin/time -v -p -o "${OUTPUT}/${chr}_kmer_gaps.fasta.time" \ + python3 "${SRCFOLDER}/2-generate_gaps/gaps_to_fasta.py" \ + "${OUTPUT}/${chr}_extracted_gaps.fa" \ + ${READSIZE} \ + "${OUTPUT}/${chr}_kmer_gaps.fasta" \ + 10 + done +} + +function align_gaps(){ + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//`; \ + bash "${SRCFOLDER}/3-align_gaps/align_gaps.sh" \ + "${OLDREF}/${chr}.${SEQ_FILE_EXT}" \ + "${OUTPUT}/${chr}_kmer_gaps.fasta" \ + ${THREAD} \ + "${OUTPUT}/${chr}_aligned_gaps" + done +} + +function extract_reads(){ + mkdir "${OUTPUT}/bedfiles" + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//`; mkdir -p "${OUTPUT}/bedfiles/${chr}/"; \ + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_regions.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_regions.sh" \ + "${OUTPUT}/${chr}_aligned_gaps.bam" \ + "${OUTPUT}/bedfiles/${chr}/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + local num_of_files=`ls -l ${OUTPUT}/bedfiles/${chr}/*.bed | wc -l` + local sed_bed_files=$(echo ${OUTPUT}\/bedfiles\/${chr}/ | sed 's/\//\\\//g') + local sed_read_scripts=$(echo ${SRCFOLDER}\/4-extract_reads\/extract_reads.bash | sed 's/\//\\\//g') + local sed_read_bam=$(echo $READ_BAM | sed 's/\//\\\//g') + cat "${SRCFOLDER}/4-extract_reads/slurm_extract_reads.sh" | \ + sed "s/{NUM_FILES}/$num_of_files/g" | \ + sed "s/{THREAD}/$THREAD/g" | \ + sed "s/{BED_FILES_FOLDER}/$sed_bed_files/g" | \ + sed "s/{EXTRACT_READS_SCRIPT}/$sed_read_scripts/g" | \ + sed "s/{READ_BAM_FILE}/$sed_read_bam/g" | \ + sed "s/{READSIZE}/$READSIZE/g" > "${OUTPUT}/${chr}_samplejob.sh" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + sbatch --wait "${OUTPUT}/${chr}_samplejob.sh" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed.time" \ + cat "${OUTPUT}/bedfiles/${chr}/"*.bed | sortBed -i - | mergeBed -i - > "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + mkdir -p "${OUTPUT}/bedfiles/${chr}/retired_bed/" + mkdir -p "${OUTPUT}/bedfiles/${chr}/constant_bed/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed.time" \ + python3 "${SRCFOLDER}/4-extract_reads/get_retired_regions.py" \ + ${READSIZE} \ + "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed" \ + $i \ + "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_reads_noprune.sh" \ + ${READ_BAM} \ + "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed" > "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + mkdir -p "${OUTPUT}/bedfiles/${chr}/reads/" + cat "${OUTPUT}/bedfiles/${chr}/"*.reads "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed" > "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed" + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_sequences_full.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_sequence.sh" \ + ${FIRST_PAIR} \ + ${SECOND_PAIR} \ + "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed" \ + "${OUTPUT}/bedfiles/${chr}/reads/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.time" \ + bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/reads_1.fastq" "${OUTPUT}/bedfiles/${chr}/reads/reads_2.fastq" | \ + samtools view -h -F4 | \ + samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.bam" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.time" \ + bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" | \ + samtools view -h -F4 | \ + samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.bam" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + bash "${SRCFOLDER}/0-align_reads.sh" \ + "${NEWREF}/$chr.${SEQ_FILE_EXT}" \ + "${OUTPUT}/bedfiles/${chr}/reads/reads" \ + "${OUTPUT}/bedfiles/${chr}/reads/paired" \ + $THREAD + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + bash "${SRCFOLDER}/0-align_singletons.sh" \ + "${NEWREF}/$chr.${SEQ_FILE_EXT}" \ + "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" \ + "${OUTPUT}/bedfiles/${chr}/reads/singletons" \ + $THREAD + done +} + +function extract_reads_with_parallel(){ + mkdir "${OUTPUT}/bedfiles" + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//`; mkdir -p "${OUTPUT}/bedfiles/${chr}/"; \ + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_regions.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_regions.sh" \ + "${OUTPUT}/${chr}_aligned_gaps.bam" \ + "${OUTPUT}/bedfiles/${chr}/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + local num_of_files=`ls -l ${OUTPUT}/bedfiles/${chr}/*.bed | wc -l` + local sed_bed_files=$(echo ${OUTPUT}\/bedfiles\/${chr}/ | sed 's/\//\\\//g') + local sed_read_scripts=$(echo ${SRCFOLDER}\/4-extract_reads\/extract_reads.bash | sed 's/\//\\\//g') + local sed_read_bam=$(echo $READ_BAM | sed 's/\//\\\//g') + cat "${SRCFOLDER}/4-extract_reads/parallel_extract_reads.sh" | \ + sed "s/{NUM_FILES}/$num_of_files/g" | \ + sed "s/{THREAD}/$THREAD/g" | \ + sed "s/{BED_FILES_FOLDER}/$sed_bed_files/g" | \ + sed "s/{EXTRACT_READS_SCRIPT}/$sed_read_scripts/g" | \ + sed "s/{READ_BAM_FILE}/$sed_read_bam/g" | \ + sed "s/{READSIZE}/$READSIZE/g" > "${OUTPUT}/${chr}_samplejob.sh" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + bash "${OUTPUT}/${chr}_samplejob.sh" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed.time" \ + cat "${OUTPUT}/bedfiles/${chr}/"*.bed | sortBed -i - | mergeBed -i - > "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + mkdir -p "${OUTPUT}/bedfiles/${chr}/retired_bed/" + mkdir -p "${OUTPUT}/bedfiles/${chr}/constant_bed/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed.time" \ + python3 "${SRCFOLDER}/4-extract_reads/get_retired_regions.py" \ + ${READSIZE} \ + "${OUTPUT}/bedfiles/${chr}/merged_${chr}.bed" \ + $i \ + "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_reads_noprune.sh" \ + ${READ_BAM} \ + "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_regions.bed" > "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + mkdir -p "${OUTPUT}/bedfiles/${chr}/reads/" + cat "${OUTPUT}/bedfiles/${chr}/"*.reads "${OUTPUT}/bedfiles/${chr}/retired_bed/retired_reads.bed" > "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed" + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/extract_sequences_full.time" \ + bash "${SRCFOLDER}/4-extract_reads/extract_sequence.sh" \ + ${FIRST_PAIR} \ + ${SECOND_PAIR} \ + "${OUTPUT}/bedfiles/${chr}/updated_and_retired_reads.bed" \ + "${OUTPUT}/bedfiles/${chr}/reads/" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.time" \ + bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/reads_1.fastq" "${OUTPUT}/bedfiles/${chr}/reads/reads_2.fastq" | \ + samtools view -h -F4 | \ + samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/paired_reads.bam" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.time" \ + bwa mem -M -t $THREAD "${NEWREF}/$chr.${SEQ_FILE_EXT}" "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" | \ + samtools view -h -F4 | \ + samtools sort -m 16g -l0 > "${OUTPUT}/bedfiles/${chr}/reads/singletons_reads.bam" + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + bash "${SRCFOLDER}/0-align_reads.sh" \ + "${NEWREF}/$chr.${SEQ_FILE_EXT}" \ + "${OUTPUT}/bedfiles/${chr}/reads/reads" \ + "${OUTPUT}/bedfiles/${chr}/reads/paired" \ + $THREAD + done + + for i in `echo ${OUTPUT}/*.chain`; do + local chr=`basename $i | sed s/.chain//` + bash "${SRCFOLDER}/0-align_singletons.sh" \ + "${NEWREF}/$chr.${SEQ_FILE_EXT}" \ + "${OUTPUT}/bedfiles/${chr}/reads/singletons.fastq" \ + "${OUTPUT}/bedfiles/${chr}/reads/singletons" \ + $THREAD + done +} + + +function main(){ + mkdir -p "${OUTPUT}/" + + local chain_existed=$(ls ${OUTPUT}/*.chain 2> /dev/null) + + if [[ ${#chain_existed} -gt 0 ]]; then + echo "Line "${LINENO}": In function "${FUNCNAME}": Use user offered chain "${chain_existed} + else + echo "Line "${LINENO}": In function "${FUNCNAME}": User does not offer chain file. Generating one." + generate_chain + fi + + generate_gaps + + align_gaps + + if [[ -z ${11} ]]; then + extract_reads_with_parallel + else + extract_reads + fi +} + + +if [[ "${#BASH_SOURCE[@]}" -eq 1 ]] +then + main "$@" +fi #constant regions: -#for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed.time" python "${SRCFOLDER}/4-extract_reads/get_constant_regions.py" ${READSIZE} $i "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed"; done -#sbatch --wrap="for i in `echo ${OUTPUT}/*.chain`; do chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_reads.bed.time" sh "${SRCFOLDER}/4-extract_reads/extract_reads_noprune.sh" ${READ_BAM} "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed" > "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_reads.bed"; done" +#for i in `echo ${OUTPUT}/*.chain`; do local chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed.time" python3 "${SRCFOLDER}/4-extract_reads/get_constant_regions.py" ${READSIZE} $i "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed"; done +#sbatch --local wrap="for i in `echo ${OUTPUT}/*.chain`; do local chr=`basename $i | sed s/.chain//`; /usr/bin/time -v -p -o "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_reads.bed.time" bash "${SRCFOLDER}/4-extract_reads/extract_reads_noprune.sh" ${READ_BAM} "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_regions.bed" > "${OUTPUT}/bedfiles/${chr}/constant_bed/constant_reads.bed"; done"