Efficient bam-to-mnd

*****************************************************

Simplified mega script for ENCODE DCC hic-pipeline.

USAGE: ./mega.sh [-g|--genome-id genome_id] [-c|--chrom-sizes <path_to_chrom.sizes>] [-v|--vcf <path_to_vcf>] [-p|--psf <path_to_psf>] [--reads-to-homologs <path_to_reads_to_homologs_file>] <path_to_merged_dedup_bam_1> ... <path_to_merged_dedup_bam_N>

USAGE: ./mega.sh -c|--chrom-sizes <path_to_chrom_sizes_file> [-g|--genome-id genome_id] [-r|--resolutions resolutions_string] [-v|--vcf <path_to_vcf>] [-C|--exclude-chr-from-diploid chr_list]  [--separate-homologs] [-p|--psf <path_to_psf>] [--reads-to-homologs <path_to_reads_to_homologs_file>] [--juicer-dir <path_to_juicer_dir>] [--phaser-dir <path_to_phaser_dir>] [-t|--threads thread_count] [-T|--threads-hic hic_thread_count] [--from-stage stage] [--to-stage stage] <path_to_merged_dedup_bam_1> ... <path_to_merged_dedup_bam_N>

DESCRIPTION:

This is a simplied mega.sh script to produce aggregate Hi-C maps and chromatic accessibility tracks from multiple experiments. The pipeline includes optional diploid steps which produce diploid versions of the contact maps and chromatin accessibility tracks.

						Shows this help.

HAPLOID PORTION

-g|--chrom-sizes    [path_to_chrom.sizes_or_genome_id]

                         Path to chrom.sizes file for the reference used when processing the individual Hi-C experiments. For several common genomes can be substituted for \"genome id\": hg18, hg19, hg38, mm9, mm10, GRCh38.

-g|--genome-id [genome_id]

                        Genome id, e.g. hg38, for some of the common references used by Juicer. Used to run the motif finder.

-c|--chrom-sizes [path_to_chrom_sizes_file]                         

                        Path to chrom.sizes file for the reference used when processing the individual Hi-C experiments.

-r|--resolutions    [string]

                        Comma-separated resolutions at which to build the hic files. Default: 2500000,1000000,500000,250000,100000,50000,25000,10000,5000,2000,1000,500,200,100,50,20,10.

DIPLOID PORTION:

-v|--vcf [path_to_vcf]

						Path to a Variant Call Format (vcf) file containing phased sequence variation data, e.g. as generated by the ENCODE DCC Hi-C variant calling & phasing pipeline.

						Path to a Variant Call Format (vcf) file containing phased sequence variation data, e.g. as generated by the ENCODE DCC Hi-C variant calling & phasing pipeline. Passing a vcf file invokes a diploid section of the script.

-reads-to-homologs [path_to_reads_to_homologs_file]

                        Path to a reads_to_homologs file generated by the phasing pipeline. If the same bams were used in the phasing pipeline this greatly reduces the compute necessary for diploid processing.

-C|--exclude-chr-from-diploid [chr_list]

						Remove specific molecules from the chromosome list (default: chrY). Note that -c and -C are incompatible in that if both are invoked, only the last listed option will be taken into account. Deafult: chrY.

--separate-homologs

						Build two separate contact maps & two separate accessibility tracks for homologs as opposed to a single diploid contact map with interleaved homologous chromosomes. This is the preferred mode for comparing contacts across homologs using the \"Observed over Control\" view for map comparison in Juicebox. Importantly, assignment of chromosomes to homologs is arbitrary. Default: not invoked.

-p|--psf [path_to_psf]

						Path to a 3D-DNA phaser psf file containing heterozygous variant and phasing information. Optional input to fast-forward some steps & save compute on processing the vcf file.

--reads-to-homologs [path_to_reads_to_homologs_file]

                        Path to a reads_to_homologs file generated by the phasing pipeline. Optional input to fast-forward some steps & save compute on diploid processing (assumes the input bams were used for phasing).

WORKFLOW CONTROL:

-t|--threads [num]

"

# defaults:

genome_id=hg38 #[talk with Muhammad, most will not work with both??]

chromSizes=hg38.chrom.sizes

resolutionsToBuildString="-r 2500000,1000000,500000,250000,100000,50000,25000,10000,5000,2000,1000,500,200,100,50,20,10"

exclude_chr="Y|chrY|MT|chrM"

separate_homolog_maps=false

mapq=1 ##lowest mapq of interest

# multithreading

threads=`parallel --number-of-cores`

tmp=`awk '/MemTotal/ {threads=int($2/1024/1024/2/6-1)}END{print threads+0}' /proc/meminfo 2>/dev/null`

tmp=$((tmp+0))

([ $tmp -gt 0 ] && [ $tmp -lt $threads ]) && threads=$tmp

threadsHic=24

[ $threadsHic -gt $threads ]&& threadsHic=$threads

threadsHic=$threads

## Create temporary directory for parallelization

tmpdir="HIC_tmp"

if [ ! -d "$tmpdir" ]; then

    mkdir "$tmpdir"

    chmod 777 "$tmpdir"

fi

export TMPDIR=${tmpdir}

# tmpdir="HIC_tmp"

# if [ ! -d "$tmpdir" ]; then

#     mkdir "$tmpdir"

#     chmod 777 "$tmpdir"

# fi

# export TMPDIR=${tmpdir}

#staging

first_stage="prep"

			exit 0

        ;;

## HAPLOID PORTION        

        -g|--chrom-sizes) OPTARG=$2 ###TODO: TALK WITH MUHAMMAD

        -g|--genome-id) OPTARG=$2

            genome_id=$OPTARG

            shift

        ;;

        -c|--chrom-sizes) OPTARG=$2

            if [ -s $OPTARG ] && [[ $OPTARG == *.chrom.sizes ]]; then

                echo "... -c|--chrom-sizes flag was triggered with $OPTARG value." >&1

                chrom_sizes=$OPTARG

            else

                echo " :( Chrom.sizes file is not found at expected location, is empty or does not have the expected extension. Exiting!" >&2

                exit 1

            fi            

            shift

        ;;

        -r|--resolutions) OPTARG=$2

            resolutionsToBuildString="-r "$OPTARG

            shift

        ;;

## DIPLOID PORTION

        -v|--vcf) OPTARG=$2

            if [ -s $OPTARG ] && [ "$OPTARG" == "*.vcf" ]; then

            if [ -s $OPTARG ] && [[ $OPTARG == *.vcf ]]; then

                echo "... -v|--vcf flag was triggered, will try to generate diploid versions of the hic file and accessibility track based on phasing data in $OPTARG." >&1

                vcf=$OPTARG

            else

                	echo " :( File is not found at expected location, is empty or does not have the expected extension. Exiting!" >&2

                	echo " :( Vcf file is not found at expected location, is empty or does not have the expected extension. Exiting!" >&2

					exit 1

            fi            

            shift

        ;;

        -C|--exclude-chr-from-diploid) OPTARG=$2

			echo "... -C|--exclude-chr-from-diploid flag was triggered, will ignore variants on $OPTARG." >&1

			exclude_chr=$OPTARG

        	shift

        ;;

        --separate-homologs)

			echo "... --separate-homologs flag was triggered, will build two separate contact maps and two separate accessibility tracks (-r.hic and -a.hic) for chromosomal homologs with identical chromosomal labels." >&1

			separate_homologs=true

		;;

        -p|--psf) OPTARG=$2

            if [ -s $OPTARG ] && [[ $OPTARG == *.psf ]]; then

                echo "... -p|--psf flag was triggered, will try to generate diploid versions of the hic file and accessibility track based on phasing data in $OPTARG." >&1

                psf=$OPTARG

            else

                	echo " :( Psf file is not found at expected location, is empty or does not have the expected extension. Exiting!" >&2

					exit 1

            fi            

            shift

        ;;

        --reads-to-homologs) OPTARG=$2

            if [ -s $OPTARG ] && [ "$OPTARG" == "*.txt" ]; then

                echo "... -r|--rth flag was triggered, will try to generate diploid versions of the hic file and accessibility track based on reads-to-homolog data in $OPTARG." >&1

            if [ -s $OPTARG ] && [[ $OPTARG == *.txt ]]; then

                echo "... --reads-to-homologs flag was triggered, will try to generate diploid versions of the hic file and accessibility track based on reads-to-homolog data in $OPTARG." >&1

                reads_to_homologs=$OPTARG

            else

                	echo " :( File is not found at expected location, is empty or does not have the expected extension. Exiting!" >&2

            fi      

            shift

        ;;

        --separate-homologs)

			echo "... --separate-homologs flag was triggered, will build two separate contact maps and two separate accessibility tracks (-r.hic and -a.hic) for chromosomal homologs with identical chromosomal labels." >&1

			separate_homologs=true

		;;

## WORKFLOW

        -t|--threads) OPTARG=$2

        	re='^[0-9]+$'

	shift

done

## TODO: give error if diploid options are invoked without a vcf file

## TODO: check that threads requested is reasonable etc.

## optional TODO: give error if diploid options are invoked without a vcf file

if [[ "${stage[$first_stage]}" -gt "${stage[$last_stage]}" ]]; then

	echo >&2 ":( Please make sure that the first stage requested is in fact an earlier stage of the pipeline to the one requested as last. Exiting!"

	exit 1

fi

[ -z $chrom_sizes ] && { echo >&2 ":( Chrom.sizes file is not optional. Please use the -c flag to point to the chrom.sizes file used when running Juicer. Exiting!"; exit 1; }

[ -z $genome_id ] && { echo >&2 ":| Warning: no genome_id is provided. Please provide a genome_id if using one of the common references to be able to run the motif finder. Ignoring motif finder!"; }

############### HANDLE DEPENDENCIES ###############

## Juicer & Phaser

bam=`echo "${@:1}"`

##TODO: check file extentions

############### MAIN #################

## 0. PREP BAM FILE

    export SHELL=$(type -p bash)

    doit () {

            samtools view -@ 2 -h reads.sorted.bam $1 | awk -F '\t' 'BEGIN{OFS="\t"}$7=="="{for(i=12;i<=NF;i++){if($i~/^ip:i:/){ip=substr($i,6)}else if ($i~/^mp:i:/){mp=substr($i,6)}else if ($i~/^MQ:i:/){mq=substr($i,6)}}; if(mq<1){next}; if(ip<mp){print}else if (ip==mp){keep($1)=$0}}' | sort -k1,1 -k2,2n -S 6G

            mapq=$2

            samtools view -@ 2 -h -q $mapq reads.sorted.bam $1 | awk -F '\t' -v mapq=$mapq '{for(i=12;i<=NF;i++){if($i~/^ip:i:/){ip=substr($i,6)}else if ($i~/^mp:i:/){mp=substr($i,6)}else if ($i~/^MQ:i:/){mq=substr($i,6)}}}(mq<mapq){next}$7=="="{if(ip>mp){next}else if (ip==mp){keep[$1]=$0}else{print 0, $3, ip, 0, 0, $3, mp, 1};next}$7<$3{next}{print 0, $3, ip, 0, 0, $7, mp, 1 > "/dev/stderr"}END{for(i in keep){n=split(keep[i],a,"\t"); for(s=12;s<=n;s++){if(a[s]~"^ip:i:"){ip=substr(a[s],6)}}; print 0, a[3], ip, 0, 0, a[3], ip, 1}}'

    }

    export -f doit

    awk '{print $1}' $chromSizes | parallel -j $threads --will-cite --joblog temp.log -k doit | sort -k1,1 -k2,2n -S 6G > tmp.bedgraph

    ##this will not work

    samtools view -@ $threads -F 1024 -O sam reads.sorted.bam | awk -v mapq=1 -f "${juicer_dir}"/scripts/sam_to_pre.awk > merged1.txt

    [ `echo "${PIPESTATUS[@]}" | tr -s ' ' + | bc` -eq 0 ] || { echo ":( Pipeline failed at sam_to_pre. Check stderr for more info. Exiting!" | tee -a /dev/stderr && exit 1; }

    ## opt1

    ## extra variable contains all (small) sequences that are not already in the chrom.sizes file. There is no use in them for generating the hic file, they are generated only for consistensy (and potentially stats)

    extra=`samtools view -H reads.sorted.bam | grep '^@SQ' | sed "s/.*SN:\([^\t]*\).*/\1/g" | awk 'FILENAME==ARGV[1]{drop[$1]=1;next}!($1 in drop)' ${chrom_sizes} - | xargs`

    #merged1.txt

    awk -v extra="$extra" '{print $1}END{print extra}' $chrom_sizes | parallel -j $threads --will-cite --joblog temp.log -k doit {} 1 >merged1.txt 2>merged1.tmp.txt

    exitval=`awk 'NR>1{if($7!=0){c=1; exit}}END{print c+0}' temp.log`

	[ $exitval -eq 0 ] || { echo ":( Pipeline failed at generating the mega mnd file. See stderr for more info. Exiting! " | tee -a /dev/stderr && exit 1; }

    rm temp.log

    sort -k2,2 -k6,6 -S 6G --parallel=${threads} merged1.tmp.txt >> merged1.txt && rm merged1.tmp.txt

    #merged30.txt

    ## MUHAMMAD: why do we have two now rather than running pre from the same one?

    awk -v extra="$extra" '{print $1}END{print extra}' $chrom_sizes | parallel -j $threads --will-cite --joblog temp.log -k doit {} 30 >merged30.txt 2>merged30.tmp.txt

    exitval=`awk 'NR>1{if($7!=0){c=1; exit}}END{print c+0}' temp.log`

	[ $exitval -eq 0 ] || { echo ":( Pipeline failed at generating the mega mnd file. See stderr for more info. Exiting! " | tee -a /dev/stderr && exit 1; }

    rm temp.log

    sort -k2,2 -k6,6 -S 6G --parallel=${threads} merged30.tmp.txt >> merged30.txt && rm merged30.tmp.txt

    ##this will not work

    samtools view -@ $threads -F 1024 -O sam reads.sorted.bam | awk -v mapq=30 -f "${juicer_dir}"/scripts/sam_to_pre.awk > merged30.txt

	[ `echo "${PIPESTATUS[@]}" | tr -s ' ' + | bc` -eq 0 ] || { echo ":( Pipeline failed at sam_to_pre. Check stderr for more info. Exiting!" | tee -a /dev/stderr && exit 1; }

    ## opt2

    # awk -v extra=$extra '{print $1}END{print extra}' $chrom_sizes | parallel -j $threads --will-cite --joblog temp.log "doit {} 1 >out.{#}.txt 2>err.{#}.txt"

    # exitval=`awk 'NR>1{if($7!=0){c=1; exit}}END{print c+0}' temp.log`

	# [ $exitval -eq 0 ] || { echo ":( Pipeline failed at building diploid contact maps. See stderr for more info. Exiting! " | tee -a /dev/stderr && exit 1; }

	# awk '{print NR}END{print NR+1}' $chrom_sizes | parallel --will-cite -k "cat out.{}.txt" >merged1.txt

    # awk '{print NR}END{print NR+1}' $chrom_sizes | parallel --will-cite -k "sort -k6,6 -S6G err.{}.txt" >>merged1.txt

    # awk '{print NR}END{print NR+1}' $chrom_sizes | parallel --will-cite "rm out.{}.txt err.{}.txt"

    # rm temp.log

    ## STATISTICS???? for now dummy files

    touch inter.txt inter_hists.m inter_30.txt inter_hists_30.m

    if [[ $threadsHic -gt 1 ]] && [[ ! -s merged1_index.txt ]]

    then

        "${juicer_dir}"/scripts/common/index_by_chr.awk merged1.txt 500000 > merged1_index.txt

        "${juicer_dir}"/scripts/index_by_chr.awk merged1.txt 500000 > merged1_index.txt

        tempdirPre="HIC_tmp" && mkdir $tempdirPre

        threadHicString="--threads $threadsHic -i merged1_index.txt -t ${tempdirPre}"

    else

        threadHicString=""

    fi

    "${juicer_dir}"/scripts/common/juicer_tools pre -n -s inter.txt -g inter_hists.m -q 1 "$resolutionsToBuildString" "$threadHicString" merged1.txt inter.hic "$genomeID"

    "${juicer_dir}"/scripts/common/juicer_tools addNorm --threads $threadsHic inter.hic

	rm -Rf "${tempdirPre}"

    ## TODO: check for failure

    if [[ $threadsHic -gt 1 ]] && [[ ! -s merged30_index.txt ]]

	then

	    "${juicer_dir}"/scripts/common/index_by_chr.awk merged30.txt 500000 > merged30_index.txt

        threadHic30String=""

	fi

    "${juicer_dir}"/scripts/common/juicer_tools pre -n -s inter_30.txt -g inter_30_hists.m -q 30 "$resolutionsToBuildString" "$threadHic30String" merged30.txt inter_30.hic "$genomeID"

	"${juicer_dir}"/scripts/common/juicer_tools addNorm --threads $threadsHic "${outputDir}"/inter_30.hic

    rm -Rf "${tempdirPre30}"

    ##MUHAMMAD:

    export IBM_JAVA_OPTIONS="-Xmx150000m -Xgcthreads1"

    export _JAVA_OPTIONS="-Xmx150000m -Xms150000m"

    ##MUHAMMAD:

    "${juicer_dir}"/scripts/juicer_tools pre -n -s inter.txt -g inter_hists.m -q 1 "$resolutionsToBuildString" "$threadHicString" merged1.txt inter.hic "$chrom_sizes"

    "${juicer_dir}"/scripts/juicer_tools addNorm --threads $threadsHic inter.hic

	rm -Rf "${tempdirPre}"

    ## TODO: check for failure

    "${juicer_dir}"/scripts/juicer_tools pre -n -s inter_30.txt -g inter_30_hists.m -q 30 "$resolutionsToBuildString" "$threadHic30String" merged30.txt inter_30.hic "$chrom_sizes"

	"${juicer_dir}"/scripts/juicer_tools addNorm --threads $threadsHic "${outputDir}"/inter_30.hic

    rm -Rf "${tempdirPre30}"

    ## TODO: check for failure

    echo ":) Done building mega hic files." >&1

    }

    export -f doit

    awk '{print $1}' $chromSizes | parallel -j $threads --will-cite --joblog temp.log -k doit | sort -k1,1 -k2,2n -S 6G > tmp.bedgraph

    awk '{print $1}' $chrom_sizes | parallel -j $threads --will-cite --joblog temp.log -k doit | sort -k1,1 -k2,2n -S 6G > tmp.bedgraph

    exitval=`awk 'NR>1{if($7!=0){c=1; exit}}END{print c+0}' temp.log`

	[ $exitval -eq 0 ] || { echo ":( Pipeline failed at building diploid contact maps. See stderr for more info. Exiting! " | tee -a /dev/stderr && exit 1; }

	rm temp.log

    bedGraphToBigWig tmp.bedgraph $chromSizes inter.bw

    bedGraphToBigWig tmp.bedgraph $chrom_sizes inter.bw

    echo ":) Done building accessibility tracks." >&1

    if [ -z $reads_to_homologs ]; then

        chr=$(awk '{str=str"|"$1}END{print substr(str,2)}' ${chromSizes})

        chr=$(awk '{str=str"|"$1}END{print substr(str,2)}' ${chrom_sizes})

        if [ -z $psf ]; then

            awk -v chr=${chr} -v output_prefix="out" -f ${phaser_dir}/phase/vcf-to-psf-and-assembly.awk ${vcf}

#     }

#     export -f doit

#     awk '{print $1}' $chromSizes | parallel -j $threads --will-cite --joblog temp.log -k doit > tmp.bedgraph

#     awk '{print $1}' $chrom_sizes | parallel -j $threads --will-cite --joblog temp.log -k doit > tmp.bedgraph

#     exitval=`awk 'NR>1{if($7!=0){c=1; exit}}END{print c+0}' temp.log`

# 	[ $exitval -eq 0 ] || { echo ":( Pipeline failed at building diploid contact maps. See stderr for more info. Exiting! " | tee -a /dev/stderr && exit 1; }