comparison libs/BWA_analysis_O_new_dependent.py @ 7:3d6680af0bec draft

planemo upload commit a4fb57231f274270afbfebd47f67df05babffa4a-dirty

6:b6281a377a18

    for_core_id=forword_seq.split(".fastq")[0]
    re_core_id=reverse_seq.split(".fastq")[0]
    for_fq=for_core_id+".fastq"
    rev_fq=re_core_id+".fastq"
    dirpath = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))#######03152016
    print "check fastq id and make them in accordance with each other...please wait..."
    os.system("python "+dirpath+"/compare_and_change_two_fastq_id.py "+for_fq+" "+rev_fq)#######03152016
    for_sai=for_core_id+".sai"
    rev_sai=re_core_id+".sai"
    sam=for_core_id+".sam"
    bam=sam.replace(".sam",".bam")

      name_list.append(line.split("\t")[2])
  a,b=Uniq(name_list)
  c=dict(zip(a,b))
  final_O=sorted(c.iteritems(), key=lambda d:d[1], reverse = True) #order from frequency high to low, but tuple while not list
  Sero_list_O=[]
  print "Final_Otype_list:"
  print final_O
  num_1=0#new inserted
  O9_wbav=0
  O310_wzx=0
  O946_wzy=0
  if len(final_O)>0: #new inserted

        O94627=x[1]        
  O_list=[]
  O_choice=""


  print "$$$Genome:",sra_name
  if len(final_O)==0:
    print "$$$No Otype, due to no hit"
  else:
    if final_O[0][1]<8:
      print "$$$No Otype, due to the hit reads number is small."
    else:
      for x in final_O:
        if x[1]>5:
          O_list.append(x[0])
      qq=1#
      for x in final_O:#
        if "sdf" in x[0] and x[1]>3:#
          qq=0#
          print "$$$",x[0],"got a hit, reads:",x[1]#
      if qq!=0:#
        print "$$$No sdf exists"#

      if "O-9,46_wbaV" in O_list and float(O9_wbaV)/float(num_1) > 0.1:
        if "O-9,46_wzy" in O_list and float(O946_wzy)/float(num_1) > 0.1:
          O_choice="O-9,46"
          print "$$$Most possilble Otype:  O-9,46"
        elif "O-9,46,27_partial_wzy" in O_list and float(O94627)/float(num_1) > 0.1:
          O_choice="O-9,46,27"
          print "$$$Most possilble Otype:  O-9,46,27"
        else:
          O_choice="O-9"
          if file_mode=="3":
            rev_fq=""
            rev_sai=""

          else:
            assembly(sra_name,O_choice,for_fq,rev_fq,for_sai,rev_sai,sam,bam,mapping_mode)
      elif ("O-3,10_wzx" in O_list) and ("O-9,46_wzy" in O_list) and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1:
        if "O-3,10_not_in_1,3,19" in O_list and float(O310_no_1319)/float(num_1) > 0.1:
          O_choice="O-3,10"
          print "$$$Most possilble Otype:  O-3,10"
        else:
          O_choice="O-1,3,19"
          print "$$$Most possilble Otype:  O-1,3,19"
      else:
        try: 
          O_choice=final_O[0][0].split("_")[0]
          if O_choice=="O-1,3,19":
            O_choice=final_O[1][0].split("_")[0]
          print "$$$Most possilble Otype: ",O_choice
        except:
          print "$$$No suitable Otype, or failure of mapping (please check the quality of raw reads)"


def assembly(sra_name,potential_choice,for_fq,rev_fq,for_sai,rev_sai,sam,bam,mapping_mode):
  database="ParaA_rfb.fasta"
  os.system(BwaPath+" index database/"+database)###01/28/2015

        O9_bigger+=1
      if O9_score<O2_score:
        O2_bigger+=1
    except:
      continue
  print "$$$Genome:",sra_name
  if O9_bigger>O2_bigger:
    print "$$$Most possible Otype is O-9"
  elif O9_bigger<O2_bigger:
    print "$$$Most possible Otype is O-2"
  else:
    print "$$$No suitable one, because can't distinct it's O-9 or O-2, but ",potential_choice," has a more possibility."
  print "O-9 number is:",O9_bigger
  print "O-2 number is:",O2_bigger

  os.system("rm "+sam+"_title.txt")###01/28/2015
  os.system("rm "+sam+"_seq.txt")###01/28/2015
  os.system("rm "+sam+".fasta")###01/28/2015
  os.system("rm "+database2+"_db.*")###01/28/2015

7:3d6680af0bec

    for_core_id=forword_seq.split(".fastq")[0]
    re_core_id=reverse_seq.split(".fastq")[0]
    for_fq=for_core_id+".fastq"
    rev_fq=re_core_id+".fastq"
    dirpath = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))#######03152016
    print("check fastq id and make them in accordance with each other...please wait...")
    os.system("python "+dirpath+"/compare_and_change_two_fastq_id.py "+for_fq+" "+rev_fq)#######03152016
    for_sai=for_core_id+".sai"
    rev_sai=re_core_id+".sai"
    sam=for_core_id+".sam"
    bam=sam.replace(".sam",".bam")

      name_list.append(line.split("\t")[2])
  a,b=Uniq(name_list)
  c=dict(zip(a,b))
  final_O=sorted(c.iteritems(), key=lambda d:d[1], reverse = True) #order from frequency high to low, but tuple while not list
  Sero_list_O=[]
  print("Final_Otype_list:")
  print(final_O)
  num_1=0#new inserted
  O9_wbav=0
  O310_wzx=0
  O946_wzy=0
  if len(final_O)>0: #new inserted

        O94627=x[1]        
  O_list=[]
  O_choice=""


  print("$$$Genome:",sra_name)
  if len(final_O)==0:
    print("$$$No Otype, due to no hit")
  else:
    if final_O[0][1]<8:
      print("$$$No Otype, due to the hit reads number is small.")
    else:
      for x in final_O:
        if x[1]>5:
          O_list.append(x[0])
      qq=1#
      for x in final_O:#
        if "sdf" in x[0] and x[1]>3:#
          qq=0#
          print("$$$",x[0],"got a hit, reads:",x[1])#
      if qq!=0:#
        print("$$$No sdf exists")#

      if "O-9,46_wbaV" in O_list and float(O9_wbaV)/float(num_1) > 0.1:
        if "O-9,46_wzy" in O_list and float(O946_wzy)/float(num_1) > 0.1:
          O_choice="O-9,46"
          print("$$$Most possilble Otype:  O-9,46")
        elif "O-9,46,27_partial_wzy" in O_list and float(O94627)/float(num_1) > 0.1:
          O_choice="O-9,46,27"
          print("$$$Most possilble Otype:  O-9,46,27")
        else:
          O_choice="O-9"
          if file_mode=="3":
            rev_fq=""
            rev_sai=""

          else:
            assembly(sra_name,O_choice,for_fq,rev_fq,for_sai,rev_sai,sam,bam,mapping_mode)
      elif ("O-3,10_wzx" in O_list) and ("O-9,46_wzy" in O_list) and float(O310_wzx)/float(num_1) > 0.1 and float(O946_wzy)/float(num_1) > 0.1:
        if "O-3,10_not_in_1,3,19" in O_list and float(O310_no_1319)/float(num_1) > 0.1:
          O_choice="O-3,10"
          print("$$$Most possilble Otype:  O-3,10")
        else:
          O_choice="O-1,3,19"
          print("$$$Most possilble Otype:  O-1,3,19")
      else:
        try: 
          O_choice=final_O[0][0].split("_")[0]
          if O_choice=="O-1,3,19":
            O_choice=final_O[1][0].split("_")[0]
          print("$$$Most possilble Otype: ",O_choice)
        except:
          print("$$$No suitable Otype, or failure of mapping (please check the quality of raw reads)")


def assembly(sra_name,potential_choice,for_fq,rev_fq,for_sai,rev_sai,sam,bam,mapping_mode):
  database="ParaA_rfb.fasta"
  os.system(BwaPath+" index database/"+database)###01/28/2015

        O9_bigger+=1
      if O9_score<O2_score:
        O2_bigger+=1
    except:
      continue
  print("$$$Genome:",sra_name)
  if O9_bigger>O2_bigger:
    print("$$$Most possible Otype is O-9")
  elif O9_bigger<O2_bigger:
    print("$$$Most possible Otype is O-2")
  else:
    print("$$$No suitable one, because can't distinct it's O-9 or O-2, but ",potential_choice," has a more possibility.")
  print("O-9 number is:",O9_bigger)
  print("O-2 number is:",O2_bigger)

  os.system("rm "+sam+"_title.txt")###01/28/2015
  os.system("rm "+sam+"_seq.txt")###01/28/2015
  os.system("rm "+sam+".fasta")###01/28/2015
  os.system("rm "+database2+"_db.*")###01/28/2015