[DEV] hicpro2juicebox now support old allValidPairs format

Author: nservant

NEWS

@@ -3,6 +3,8 @@ CHANGES IN VERSION 2.10.0
 
 NEW FEATURES
 
+   o hicpro2juicebox.sh utility now supports alof HiC-Pro format (< 2.7.5)
+
 SIGNIFICANT USER-VISIBLE CHANGES
 
    o udpate R scripts to be compatible with the lastest ggplot2 version (>2.2.1)
@@ -11,6 +13,8 @@ SIGNIFICANT USER-VISIBLE CHANGES
 
 BUG FIXES
 
+   o Fix bug to avoid floating values in valid pair positions 
+
    o Fix bug in order of samtools sort parameter in bowtie_combine.sh
 
    o Fix bug in output option of makeViewpoints

bin/utils/extract_snps.py

@@ -28,6 +28,7 @@ def usage():
     print "-a/--alt <sample name for alternative allele>"
     print "[-r/--ref] <sample name for reference allele. Default is the same as in the initial file>"
     print "[-f/--filt] <Filtering level. 0 - no filtering. 1 - based on FI information. 2 - based on FILTER information. Default is 2>"
+    print "[-x/--exclude] <Exclude potential contaminants from the list of SNPs which are equal to ALT genotype, i.e REF is expected to be contaminated>"
     print "[-v/--verbose] <Verbose>"
     print "[-h/--help] <Help>"
     return
@@ -37,11 +38,12 @@ def get_args():
     try:
         opts, args = getopt.getopt(
             sys.argv[1:],
-            "i:r:a:f:vh",
+            "i:r:a:f:x:vh",
             ["vcf=",
              "alt=",
              "ref=",
              "filt=", 
+             "excl=",
              "verbose", "help"])
     except getopt.GetoptError:
         usage()
@@ -50,7 +52,7 @@ def get_args():
 
 ## 0/0, 1/1, A, T, 1
 ## 1/1, 2/2, A, T, C, 1
-def get_filter_snp_gt(gref, galt, ref, alt):
+def get_filter_snp_gt(gref, galt, ref, alt, conta):
 
     #print gref
     #print galt
@@ -76,14 +78,21 @@ def get_filter_snp_gt(gref, galt, ref, alt):
         ## Non informative SNPs
         if ref_snp == alt_snp:
             return -3
-        else:
-            ## check alt and ref alleles
-            alleles = []
-            alleles.append(ref)
-            for a in alt.split(','):
-                alleles.append(a)
+        ## Remove SNPs that are different from the reference and therefore likely to be wrongly assigned to the alternative
+        elif len(conta)>0:
+            for i in range(len(conta)):
+                conta_geno = re.split('/|\|', conta[i])
+                ## if conta = alt remove snps
+                if conta_geno[0] == alt_geno[0] or conta_geno[1] == alt_geno[1]:
+                    return -4
+        
+        ## check alt and ref alleles
+        alleles = []
+        alleles.append(ref)
+        for a in alt.split(','):
+            alleles.append(a)
 
-            return [alleles[int(ref_snp)], alleles[int(alt_snp)]]
+        return [alleles[int(ref_snp)], alleles[int(alt_snp)]]
 
 
 if __name__ == "__main__":
@@ -93,6 +102,7 @@ def get_filter_snp_gt(gref, galt, ref, alt):
     vcfFile = None
     refSample = None
     altSample = None
+    exclusion = None
     filt_qual=2
     verbose = False
 
@@ -112,6 +122,8 @@ def get_filter_snp_gt(gref, galt, ref, alt):
             altSample = arg
         elif opt in ("-f", "--filt"):
             filt_qual = int(arg)
+        elif opt in ("-x", "--exclude"):
+            exclusion = arg
         elif opt in ("-v", "--verbose"):
                 verbose = True
         else:
@@ -133,16 +145,20 @@ def get_filter_snp_gt(gref, galt, ref, alt):
     samples = []
     altidx = -1
     refidx = -1
-
+    contaidx=[]
+    
     var_counter = 0
     snp_counter = 0
     hetero_counter = 0
     badqual_counter = 0
     undefined_counter = 0
     nonspe_counter = 0
+    conta_counter = 0
 
     for line in vcf_handle:
         line = line.rstrip()
+        #print >> sys.stderr, line
+
         ## for now we don't care about the header 
         if line.startswith('##'):
             if refSample is not None and line.startswith("##reference="):
@@ -158,7 +174,29 @@ def get_filter_snp_gt(gref, galt, ref, alt):
                     refidx = i
                 elif samples[i] == altSample:
                     altidx = i
-                    
+                elif exclusion is not None:
+                    ## conta idx
+                    exs=exclusion.split(",")
+                    for i in range(len(exs)):
+                        ct = exs[i]
+                        if samples[i] == ct:
+                            contaidx.append(i)
+                            if verbose:
+                                print >> sys.stderr, "## Potential Contaminant(s) = " + ct                           
+
+
+            ## Check if Bl6 is in the conta list
+            if exclusion is not None:
+                exs=exclusion.split(",")
+                for i in range(len(exs)):
+                    ct = exs[i]
+                    if ct == "REF":
+                        contaidx.append(-1)
+                        if verbose:
+                            print >> sys.stderr, "## Potential Contaminant(s) = REF"
+
+
+
             ## Check input parameters
             if refSample != None and refidx == -1:
                 print   >> sys.stderr, "Error : REF sample not found"
@@ -185,7 +223,10 @@ def get_filter_snp_gt(gref, galt, ref, alt):
 
             fields = line.split('\t',9)
             var_counter+=1
-            
+
+            ## init list of contaminant
+            contg=[]
+
             ## check chromosomes name
             if re.compile('^chr').match(fields[0]):
                 chrom=fields[0]
@@ -214,24 +255,37 @@ def get_filter_snp_gt(gref, galt, ref, alt):
                 else:
                     refg = ["0/0"]
                     reffi =  "1"
+                    
+                if len(contaidx) > 0:
+                    for  i in range(len(contaidx)):
+                        if contaidx[i] == -1:
+                            contg.append("0/0")
+                        else:
+                            cg=genotypes[contaidx[i]].split(':')
+                            cfi = cg[len(cg)-1]
+                            if filt_qual != 1 or filt_qual == 1 and cfi == str(1):
+                                contg.append(cg[0])
 
                 ## Filter on FI field
                 if filt_qual != 1 or (filt_qual == 1 and reffi == str(1) and altfi == str(1)):
                     #print "---------"
                     #print refg
                     #print altg
                     #print fields
-                    geno = get_filter_snp_gt(refg[0], altg[0], fields[3], fields[4])
+                    geno = get_filter_snp_gt(refg[0], altg[0], fields[3], fields[4], contg)
+
                     if geno == -1:
                         undefined_counter += 1 
                     elif geno == -2:
                         hetero_counter += 1
                     elif geno == -3:
                         nonspe_counter += 1
+                    elif geno == -4:
+                        conta_counter += 1
                     else:
                         snp_counter += 1
                         #altg[0]="1/1"
-                        #print chrom + "\t" + fields[1] + "\t" + fields[2] + "\t" + geno[0] + "\t" + geno[1] + "\t" + fields[5] + "\t" + fields[6] + "\t" + fields[7] + "\t" + fields[8] + "\t" + ":".join(altg)
+                        ##print chrom + "\t" + fields[1] + "\t" + fields[2] + "\t" + geno[0] + "\t" + geno[1] + "\t" + fields[5] + "\t" + fields[6] + "\t" + fields[7] + "\t" + fields[8] + "\t" + ":".join(altg)
                         print chrom + "\t" + fields[1] + "\t" + fields[2] + "\t" + geno[0] + "\t" + geno[1] + "\t" + fields[5] + "\t" + fields[6] + "\t" + fields[7] + "\t" + "GT" + "\t" + "0/1"
 
                 else:
@@ -252,6 +306,7 @@ def get_filter_snp_gt(gref, galt, ref, alt):
          print >> sys.stderr, "## Number of heterozygous SNPs =", hetero_counter
          print >> sys.stderr, "## Number of undefined genotype SNPs =", undefined_counter
          print >> sys.stderr, "## Number of bad quality SNPs =", badqual_counter
+         print >> sys.stderr, "## Number of potential contaminant SNPs =", conta_counter
 
 
     vcf_handle.close()

bin/utils/hicpro2juicebox.sh

@@ -105,6 +105,23 @@ if [[ ! -e $JUICEBOXJAR ]]; then
     exit 1
 fi
 
+## Deal with old format
+nbfields=$(head -1 $VALIDPAIRS | awk '{print NF}')
+
+if [[ $nbfields == "12" ]]; then
+    echo -e "HiC-Pro format > 2.7.5 detected ..."
+
+elif [[ $nbfields == "8" ]]; then
+    echo -e "HiC-Pro format < 2.7.6 detected ..."
+    echo -e "Adjusting AllValidPairs format ..."
+    awk '{OFS="\t"; print $0,0,1,42,42}' $VALIDPAIRS > ${TEMP}/$$_format_AllValidPairs
+    VALIDPAIRS=${TEMP}/$$_format_AllValidPairs
+
+else
+    echo -e "Error : unknown format - $nbfields detected, whereas 8 (< v2.7.6) or 12 (> v2.7.5) fields are expected !"
+    exit 1
+fi
+
 echo "Generating Juicebox input files ..."
 
 if [[ ! -z $RESFRAG ]]; then
@@ -117,7 +134,7 @@ if [[ ! -z $RESFRAG ]]; then
 
     ## The “pre” command needs the contact map to be sorted by chromosome and grouped so that all reads for one chromosome (let’s say, chr1) appear in the same column.
     ## Also, chromosomes should not have the ‘chr” substring and the strand is coded as 0 for positive and anything else for negative (in practice, 1).
-    awk '{$4=$4!="+"; $7=$7!="+"; n1=split($9, frag1, "_"); n2=split($10, frag2, "_"); } $2<=$5{print $1, $4, $2, $3, frag1[n1], $7, $5, $6, frag2[n2], $11, $12 }$5<$2{ print $1, $7, $5, $6, frag2[3], $4, $2, $3, frag1[3], $12, $11}' $VALIDPAIRS | sort -k3,3d  -k7,7d -S 90 > ${TEMP}/$$_allValidPairs.pre_juicebox_sorted
+    awk '{$4=$4!="+"; $7=$7!="+"; n1=split($9, frag1, "_"); n2=split($10, frag2, "_"); } $2<=$5{print $1, $4, $2, $3, frag1[n1], $7, $5, $6, frag2[n2], $11, $12 }$5<$2{ print $1, $7, $5, $6, frag2[n2], $4, $2, $3, frag1[n1], $12, $11}' $VALIDPAIRS | sort -k3,3d  -k7,7d -S 90 > ${TEMP}/$$_allValidPairs.pre_juicebox_sorted
 else
     awk '{$4=$4!="+"; $7=$7!="+"} $2<=$5{print $1, $4, $2, $3, 0, $7, $5, $6, 1, $11, $12 }$5<$2{ print $1, $7, $5, $6, 0, $4, $2, $3, 1, $12, $11 }' $VALIDPAIRS | sort -k3,3d  -k7,7d -S 90 > ${TEMP}/$$_allValidPairs.pre_juicebox_sorted
 fi

scripts/mapped_2hic_fragments.py

@@ -136,7 +136,7 @@ def get_read_pos(read):
     read : list
         list of aligned reads
     """
-    pos = read.pos + read.alen/2
+    pos = read.pos + int(read.alen/2)
 
     return pos