nfrDynAna2Matrix

#!/bin/bash
#PBS -l nodes=1:ppn=4

REPLICATES=2
MIN_N_CLASS=500
GENOME="mm9"
SAMPLE_COUNT=4

#### usage ####
usage() {
	echo Program: "nfrDynAna2Matrix (compute matrix containing activity features corresponding to input coordinates in BED format)"
	echo Author: BRIC, University of Copenhagen, Denmark
	echo Version: 1.0
	echo Contact: pundhir@binf.ku.dk
	echo "Usage: nfrDynAna2Matrix -i <file> -j <file> -o <dir> [OPTIONS]"
	echo "Options:"
	echo " -i <file>   [input genomic coordinates in BED format (can be stdin)]"
    echo "             [output from nfrDynAna or multiIntersectBed.sh script]"
    echo "             [chr start end stages_count stages (..)]"
    echo " -j <file>   [configuration file]"
    echo "             [<id> <bam file> (id can be: nfrdip, tpm, expr, cobound, tss)]"
    echo " -o <dir>    [output directory]"
    echo "[OPTIONS]"
    echo " -r <int>    [number of replicates. used to compute mean (default: 2)]"
    echo " -n <int>    [minimum number of elements within each nfr dynamic class (default: 500)]"
    echo " -g <string> [genome (default: mm9)]"
    echo " -a          [compute raw read counts instead (default: TPM)]"
    echo "[OPTIONS: chip2geneExprCor]"
    echo " -x <int>    [number of samples for which to compute gene expression (default: 4)]"
    echo " -c          [associate enhancers to closest gene instead to collapsed gene(s)]"
	echo " -h          [help]"
    echo "[NOTE]"
    echo " In case of error: Could not open input BAM files, refer to:"
    echo " https://groups.google.com/forum/#!msg/bedtools-discuss/_LNuoRWHn50/14MaqyzyzXsJ"
    echo " i.e. make sure header of all input bam files is identical"
	echo
	exit 0
}

#### parse options ####
while getopts i:o:j:r:n:g:ax:ch ARG; do
    case "$ARG" in
    i) BED_FILE=$OPTARG;;
    o) OUTDIR=$OPTARG;;
    j) CONFIG_FILE=$OPTARG;;
    r) REPLICATES=$OPTARG;;
    n) MIN_N_CLASS=$OPTARG;;
    g) GENOME=$OPTARG;;
    a) RAW_COUNT=1;;
    x) SAMPLE_COUNT=$OPTARG;;
    c) CLOSEST_GENE=1;;
    h) HELP=1;;
    esac
done

## usage, if necessary file and directories are given/exist
if [ -z "$BED_FILE" -o ! -f "$CONFIG_FILE" -o "$HELP" ]; then
	usage
fi

## populating files based on input genome
if [ "$GENOME" == "mm9" ]; then
    GENOME_FILE="/home/pundhir/project/genome_annotations/mouse.mm9.genome"
elif [ "$GENOME" == "hg19" ]; then
    GENOME_FILE="/home/pundhir/project/genome_annotations/human.hg19.genome"
elif [ "$GENOME" == "rn5" ]; then
    GENOME_FILE="/home/pundhir/project/genome_annotations/rat.rn5.genome"
else
    echo "Presently the program only support analysis for mm9, hg19 or rn5"
    echo
    usage
fi

## create output directory, if do not exist
if [ ! -d "$OUTDIR" ]; then
    mkdir -p $OUTDIR
fi

## create temporary BED file if input is from stdin
if [ "$BED_FILE" == "stdin" ]; then
    TMP=$(cat /dev/urandom | tr -dc 'a-zA-Z0-9' | fold -w 32 | head -n 1)
    while read LINE; do
        echo ${LINE}
    done | perl -ane '$line=""; foreach(@F) { $line.="$_\t"; } $line=~s/\t$//g; print "$line\n";' > $OUTDIR/$TMP
    BED_FILE=$OUTDIR/$TMP
fi

## filter input BED file based on minimum number of elements within each nfr dynamic class
for CLASS in $(cat $BED_FILE | cut -f 5 | sort | uniq -c | sed -E 's/^\s+//g' | sort -k 1rn,1 | perl -ane 'if($F[0]>'$MIN_N_CLASS') { print "$F[1]\n"; }'); do perl -ane 'if($F[4]=~/^'$CLASS'$/) { print $_; }' $BED_FILE; done | sort -k 1,1 -k 2n,2 -k 3n,3 > $OUTDIR/NFR_DYNAMICS_SIG
BED_FILE="NFR_DYNAMICS_SIG"

## read configuration file
BAM_FILE_NFRDIP=$(cat $CONFIG_FILE | perl -ane '
    if($_=~/^nfrdip/) {
        $file.="$F[1],";
    } END {
        $file=~s/\,$//g;
        print "$file\n";
    }'
)
#echo $BAM_FILE_NFRDIP

BAM_FILE_TPM=$(cat $CONFIG_FILE | perl -ane '
    if($_=~/^tpm/) {
        $file.="$F[1],";
    } END {
        $file=~s/\,$//g;
        print "$file\n";
    }'
)
#echo $BAM_FILE_TPM

BED_FILE_EXPR=$(cat $CONFIG_FILE | perl -ane '
    if($_=~/^expr/) {
        $file.="$F[1],";
    } END {
        $file=~s/\,$//g;
        print "$file\n";
    }'
)
#echo $BED_FILE_EXPR

BED_FILE_COBOUND=$(cat $CONFIG_FILE | perl -ane '
    if($_=~/^cobound/) {
        $file.="$F[1],";
    } END {
        $file=~s/\,$//g;
        print "$file\n";
    }'
)
#echo $BED_FILE_COBOUND

BED_FILE_TSS=$(cat $CONFIG_FILE | perl -ane '
    if($_=~/^tss/) {
        $file.="$F[1],";
    } END {
        $file=~s/\,$//g;
        print "$file\n";
    }'
)
#echo $BED_FILE_TSS

## check, if index of BAM files exist
BAM_FILE=$(echo "$BAM_FILE_NFRDIP,$BAM_FILE_TPM" | perl -ane '$_=~s/^\,//g; $_=~s/\,$//g; print $_;')
oIFS=$IFS
IFS=","
BAMFILES=($BAM_FILE)
BAMFILES_COUNT=${#BAMFILES[@]}
IFS=$oIFS

for(( i=0; i<$BAMFILES_COUNT; i++ )); do
    ## create index of input BAM file, if does not exist
    if [ ! -f "${BAMFILES[$i]}.bai" ]; then
        samtools index ${BAMFILES[$i]}
    fi
done

<<"COMMENT"
COMMENT
## compute nfr dip
NCOL=$(head $OUTDIR/$BED_FILE -n 1 | perl -ane '$ncol=scalar(@F)+1; print "$ncol";')
if [ ! -z "$BAM_FILE_NFRDIP" -a ! -s "$OUTDIR/$BED_FILE.nfrdip" ]; then
    bed2window -i $OUTDIR/$BED_FILE -w 50 -s | bed2nfrdipV2 -i stdin -j $BAM_FILE_NFRDIP -t | rep2mean -i stdin -s $NCOL -n $REPLICATES > $OUTDIR/$BED_FILE.nfrdip
else
    touch $OUTDIR/$BED_FILE.nfrdip
fi

## compute TPM for TF binding
NCOL=$(head $OUTDIR/$BED_FILE -n 1 | perl -ane '$ncol=scalar(@F)+1; print "$ncol";')
if [ ! -z "$BAM_FILE_TPM" -a ! -s "$OUTDIR/$BED_FILE.tpm" ]; then
    if [ -z "$RAW_COUNT" ]; then
        bed2expr -i $OUTDIR/$BED_FILE -j $BAM_FILE_TPM -m -d -g $GENOME | rep2mean -i stdin -s $NCOL -n $REPLICATES > $OUTDIR/$BED_FILE.tpm
        #cat $OUTDIR/$BED_FILE | bed2expr -i stdin -j $BAM_FILE_TPM -d -g $GENOME | rep2mean -i stdin -s $NCOL -n $REPLICATES > $OUTDIR/$BED_FILE.tpm
    else
        bed2expr -i $OUTDIR/$BED_FILE -j $BAM_FILE_TPM -d -g $GENOME | rep2mean -i stdin -s $NCOL -n $REPLICATES > $OUTDIR/$BED_FILE.tpm
    fi
else
    touch $OUTDIR/$BED_FILE.tpm
fi

## compute proximal gene expression
if [ ! -z "$BED_FILE_EXPR" ]; then
    if [ ! -z "$CLOSEST_GENE" ]; then
        chip2geneExprCor -i $OUTDIR/$BED_FILE -j $BED_FILE_EXPR -k $SAMPLE_COUNT -n -c > $OUTDIR/$BED_FILE.expr
    else
        chip2geneExprCor -i $OUTDIR/$BED_FILE -j $BED_FILE_EXPR -k $SAMPLE_COUNT -n > $OUTDIR/$BED_FILE.expr
    fi
else
    touch $OUTDIR/$BED_FILE.expr
fi

## compute number of proximal genes
if [ ! -z "$BED_FILE_EXPR" ]; then
    #intersectBed -a $OUTDIR/$BED_FILE -b $BED_FILE_EXPR -c > $OUTDIR/$BED_FILE.count
    paste $OUTDIR/$BED_FILE <(less $OUTDIR/$BED_FILE.expr | cut -f 5 | perl -an -F',' -e 'print scalar(@F)."\n";') > $OUTDIR/$BED_FILE.count
else
    touch $OUTDIR/$BED_FILE.count
fi

## compute distance to proximal gene
if [ ! -z "$BED_FILE_TSS" ]; then
        closestBed -a $OUTDIR/$BED_FILE -b <(sortBed -i $BED_FILE_TSS) -t first -d | perl -ane 'print $F[scalar(@F)-1]."\n";' > $OUTDIR/$BED_FILE.dist
else
    if [ ! -z "$BED_FILE_EXPR" ]; then
        closestBed -a $OUTDIR/$BED_FILE -b <(sortBed -i $BED_FILE_EXPR) -t first -d | perl -ane 'print $F[scalar(@F)-1]."\n";' > $OUTDIR/$BED_FILE.dist
    else
        touch $OUTDIR/$BED_FILE.dist
    fi
fi

## compute fraction of nfr (co)-bound by TFs
NCOL=$(head $OUTDIR/$BED_FILE -n 1 | perl -ane '$ncol=scalar(@F)+1; print "$ncol";')
if [ ! -z "$BED_FILE_COBOUND" ]; then
    MERGE_COL=$(echo $NCOL | perl -ane '$val=(); foreach($i=4; $i<'$NCOL'; $i++) { $val.="$i,"; } $val=~s/\,$//g; print "$val";')
    MERGE_COL=$(perl -e '$val='$NCOL'+4; chomp($_); print "'$MERGE_COL',$val";')
    intersectBed -a $OUTDIR/$BED_FILE -b $BED_FILE_COBOUND -wao | sortBed -i stdin | mergeBed -i stdin -c $MERGE_COL -o distinct | perl -ane '$col=scalar(@F)-1; if($F[$col]=~/pu1/ && $F[$col]=~/cebpa/) { $F[$col]="pu1,cebpa"; } elsif($F[$col]=~/pu1/) { $F[$col]="pu1"; } elsif($F[$col]=~/cebpa/) { $F[$col]="cebpa"; } $line=(); foreach(@F) { $line.="$_\t"; } $line=~s/\t$//g; print "$line\n";' > $OUTDIR/$BED_FILE.cobound
    ## R-code is doing this now
    #intersectBed -a $OUTDIR/$BED_FILE -b $BED_FILE_COBOUND -wao | sortBed -i stdin | mergeBed -i stdin -c 4,5,6,7,8,9,14 -o distinct | cut -f 5,10 | perl -ane 'if($F[1]=~/pu1/ && $F[1]=~/cebpa/) { $F[1]="pu1,cebpa"; } elsif($F[1]=~/pu1/) { $F[1]="pu1"; } elsif($F[1]=~/cebpa/) { $F[1]="cebpa"; } print "$F[0]\t$F[1]\n";' | sort | uniq -c | sed -E 's/^\s+//g' | sort -k 3,3 | perl -ane '$count{$F[1]}+=$F[0]; push(@data, $_); END { foreach(@data) { @T=split(/\s+/,$_); $fraction=sprintf("%0.6f", $_/$count{$T[1]}); chomp($_); print "$_\t$fraction\n"; } }' > $OUTDIR/$BED_FILE.cobound.summary
else
    touch $OUTDIR/$BED_FILE.cobound
fi

## join all files into one
NCOL_nfrdip=$(head $OUTDIR/$BED_FILE.nfrdip -n 1 | perl -ane '$ncol=scalar(@F); print "$ncol";')
NCOL_tpm=$(head $OUTDIR/$BED_FILE.tpm -n 1 | perl -ane '$ncol=scalar(@F); print "$ncol";')
NCOL_expr=$(head $OUTDIR/$BED_FILE.expr -n 1 | perl -ane '$ncol=scalar(@F); print "$ncol";')
NCOL_count=$(head $OUTDIR/$BED_FILE.count -n 1 | perl -ane '$ncol=scalar(@F); print "$ncol";')
NCOL_cobound=$(head $OUTDIR/$BED_FILE.cobound -n 1 | perl -ane '$ncol=scalar(@F); print "$ncol";')

paste $OUTDIR/$BED_FILE <(cut -f $NCOL-$NCOL_nfrdip $OUTDIR/$BED_FILE.nfrdip) <(cut -f $NCOL-$NCOL_tpm $OUTDIR/$BED_FILE.tpm) <(cut -f 5-$NCOL_expr $OUTDIR/$BED_FILE.expr) <(cut -f $NCOL-$NCOL_count $OUTDIR/$BED_FILE.count) <(cut -f $NCOL-$NCOL_cobound $OUTDIR/$BED_FILE.cobound) <(cat $OUTDIR/$BED_FILE.dist) | tabEdit -i stdin -r > $OUTDIR/$BED_FILE.stat
#echo -e "$NCOL\t$NCOL_nfrdip\t$NCOL_tpm\t$NCOL_expr\t$NCOL_count\t$NCOL_cobound"

## remove temporary file
if [ ! -z "$TMP" ]; then
    rm $OUTDIR/$TMP
fi