From 469ad5b6807352f22f20cba6a2dc3c3d263cdd42 Mon Sep 17 00:00:00 2001 From: yangyxt Date: Thu, 19 Nov 2020 18:11:05 +0800 Subject: [PATCH 1/5] add new way to call extract reads --- 4-extract_reads/extract_reads.sh | 12 +- 4-extract_reads/parallel_extract_reads.sh | 22 ++ 4-extract_reads/slurm_extract_reads.sh | 7 +- run_pipeline.sh | 305 ++++++++++++++++++++-- 4 files changed, 322 insertions(+), 24 deletions(-) create mode 100644 4-extract_reads/parallel_extract_reads.sh 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..5692bc6 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,290 @@ READ_BAM=$6 FIRST_PAIR=$7 SECOND_PAIR=$8 OUTPUT=$9 -THREAD=$10 +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 + +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(){ + 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 + generate_chain + fi + + generate_gaps + + align_gaps + + if [[ -z ${11} ]]; then + extract_reads_with_parallel + else + extract_reads + 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" From b654d9fd7783b83074711c69919193cb5cb2a062 Mon Sep 17 00:00:00 2001 From: yangyxt Date: Thu, 19 Nov 2020 20:11:18 +0800 Subject: [PATCH 2/5] add execution block --- run_pipeline.sh | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/run_pipeline.sh b/run_pipeline.sh index 5692bc6..757c7b0 100644 --- a/run_pipeline.sh +++ b/run_pipeline.sh @@ -294,6 +294,12 @@ function main(){ fi } + +if [[ "${#BASH_SOURCE[@]}" -eq 1 ]] +then + main "$@" +fi + #constant regions: #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" From a3142a622437e3587974a1a703bfeea9a9c96e57 Mon Sep 17 00:00:00 2001 From: yangyxt Date: Thu, 19 Nov 2020 20:27:31 +0800 Subject: [PATCH 3/5] move mkdir to main func --- run_pipeline.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/run_pipeline.sh b/run_pipeline.sh index 757c7b0..1d8660f 100644 --- a/run_pipeline.sh +++ b/run_pipeline.sh @@ -13,7 +13,6 @@ SECOND_PAIR=$8 OUTPUT=$9 THREAD=${10} -mkdir -p "${OUTPUT}/" function generate_chain(){ bash "${SRCFOLDER}/1-generate_chain/chain_install.sh" "${OUTPUT}/" @@ -275,6 +274,8 @@ function extract_reads_with_parallel(){ function main(){ + mkdir -p "${OUTPUT}/" + local chain_existed=$(ls ${OUTPUT}/*.chain 2> /dev/null) if [[ ${#chain_existed} -gt 0 ]]; then From baf4ac990bb16248a821dac1bbf68bd65fe08f82 Mon Sep 17 00:00:00 2001 From: yangyxt Date: Sun, 29 Nov 2020 22:16:54 +0800 Subject: [PATCH 4/5] add logging --- run_pipeline.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/run_pipeline.sh b/run_pipeline.sh index 1d8660f..42e97dd 100644 --- a/run_pipeline.sh +++ b/run_pipeline.sh @@ -281,6 +281,7 @@ function main(){ 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 From f06b09b71ffd45c741fab662b0dc218910295dd2 Mon Sep 17 00:00:00 2001 From: yangyxt Date: Mon, 30 Nov 2020 10:10:42 +0800 Subject: [PATCH 5/5] open file with "r" mode, otherwise will raise NotADirectory error --- 2-generate_gaps/extract_gaps.py | 2 +- run_pipeline.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/run_pipeline.sh b/run_pipeline.sh index 42e97dd..8d9b650 100644 --- a/run_pipeline.sh +++ b/run_pipeline.sh @@ -15,7 +15,7 @@ THREAD=${10} function generate_chain(){ - bash "${SRCFOLDER}/1-generate_chain/chain_install.sh" "${OUTPUT}/" + bash "${SRCFOLDER}/1-generate_chain/chain_install.sh" "${OUTPUT}" bash "${SRCFOLDER}/1-generate_chain/sample_run.sh" \ "${OLDREF}" \ "${NEWREF}" \