Python MafIndex.get_spliced方法代码示例

# 需要导入模块: from Bio.AlignIO.MafIO import MafIndex [as 别名]
# 或者: from Bio.AlignIO.MafIO.MafIndex import get_spliced [as 别名]
    class TestSpliceGoodMAF(unittest.TestCase):
        """Test in silico splicing on a correctly-formatted MAF"""

        def setUp(self):
            self.idx = MafIndex("MAF/ucsc_mm9_chr10_big.mafindex",
                                "MAF/ucsc_mm9_chr10_big.maf", "mm9.chr10")
            self.assertEqual(len(self.idx), 983)

        def test_invalid_strand(self):
            self.assertRaises(ValueError,
                              self.idx.get_spliced,
                              (0, 1000), (500, 1500), ".")

        def test_no_alignment(self):
            result = self.idx.get_spliced((0, 1000), (500, 1500), 1)

            self.assertEqual(len(result), 1)
            self.assertEqual(len(result[0].seq), 1000)
            self.assertEqual(str(result[0].seq), "N" * 1000)

        def test_correct_retrieval_1(self):
            """Correct retrieval of Cnksr3 in mouse.

            This is the real thing. We're pulling the spliced alignment of
            an actual gene (Cnksr3) in mouse. It should perfectly match the
            spliced transcript pulled independently from UCSC.
            """
            result = self.idx.get_spliced((3134303, 3185733, 3192055, 3193589,
                                           3203538, 3206102, 3208126, 3211424,
                                           3211872, 3217393, 3219697, 3220356,
                                           3225954),
                                          (3134909, 3185897, 3192258, 3193677,
                                           3203580, 3206222, 3208186, 3211493,
                                           3212019, 3217518, 3219906, 3220446,
                                           3227479), 1)

            cnksr3 = str(SeqIO.read("MAF/cnksr3.fa", "fasta").seq).upper()
            mm9_seq = "".join([str(x.seq) for x in result
                               if x.id.startswith("mm9")]).replace("-", "")

            self.assertEqual(mm9_seq, cnksr3)

# 需要导入模块: from Bio.AlignIO.MafIO import MafIndex [as 别名]
# 或者: from Bio.AlignIO.MafIO.MafIndex import get_spliced [as 别名]
def main():
  parser = argparse.ArgumentParser(description='report pair-wise PID in intervals', 
                                   prog='maf2pid')
  parser.add_argument('--splst', dest='splst', help='species list')
  parser.add_argument('--bed', dest='interval', help='bed intervals')
  parser.add_argument('--output', dest='output', help='output file')
  if len(sys.argv) < 4:
    parser.print_help()
    sys.exit(1)

  args = parser.parse_args()

  intervals_by_chrom = parse_bed(args.interval)
  [species, pairwise] = make_pair_names(args.splst)

  fout = open(args.output, 'w')
  header = "#chrom\tstart\tend\tpme"
  for i in pairwise :
    header += "\t"+i
  fout.write( header + "\n")
  
  for chrom in intervals_by_chrom.keys() :
    sqlite_file = chrom + ".mafindex"
    maf_file = chrom + ".7sp.maf"
    target_seqname = "hg19." + chrom
    idx = MafIndex(sqlite_file, maf_file, target_seqname)
    # MafIndex.get_spliced(starts, ends, strand='+1')
    starts = intervals_by_chrom[chrom][0]
    ends = intervals_by_chrom[chrom][1]
    names = intervals_by_chrom[chrom][2]
    for k in range(len(starts)) :
      result = idx.get_spliced([starts[k]], [ends[k]], strand='+1')
      recsp = [rec.id.split(".")[0] for rec in result]
      PIDscores = {}
      for i in range(len(recsp)-1) :
        for j in range(i+1, len(recsp)) :
          pw = recsp[i]+ "_"+ recsp[j]
          if pw not in pairwise :
            pw = recsp[j]+ "_"+ recsp[i]
          if (pw in pairwise) and (not PIDscores.has_key(pw)): # exclude paralogs
            pid = get_PID(str(result[i].seq), str(result[j].seq))
            PIDscores[pw] = pid
      for pw in pairwise :
        if pw not in PIDscores :
          PIDscores[pw] = -1
      fout.write(chrom + "\t" + str(starts[k]) + "\t" + str(ends[k]) + "\tpme_" + names[k])
      for pw in pairwise :
        fout.write("\t"+ "{:.5f}".format(PIDscores[pw]))
      fout.write("\n")
  fout.close()
  return 

# 需要导入模块: from Bio.AlignIO.MafIO import MafIndex [as 别名]
# 或者: from Bio.AlignIO.MafIO.MafIndex import get_spliced [as 别名]
                exon_s.append(g.initial)
                exon_e.append(g.final)
        else:
            exon_s = [rg.initial]
            exon_e = [rg.final]

    else: 
        exon_s = [rg.initial]
        exon_e = [rg.final]

    if rg.orientation == "+":
        strand = "+1"
    else:
        strand = "-1"     
    multiple_alignment = idx.get_spliced(exon_s,
                                         exon_e,
                                         strand = strand)
    
    for j, seqrec in enumerate(multiple_alignment):
        if seqrec.id.startswith(args.organism):
            if j == 0: break
            else:
                multiple_alignment._records[0], multiple_alignment._records[j] = multiple_alignment._records[j], multiple_alignment._records[0]
    
    #print(dir(multiple_alignment))
    #sys.exit(0)

    seqs = []
    for seqrec in multiple_alignment:
        #print(dir(seqrec))
        try:

# 需要导入模块: from Bio.AlignIO.MafIO import MafIndex [as 别名]
# 或者: from Bio.AlignIO.MafIO.MafIndex import get_spliced [as 别名]

#.........这里部分代码省略.........
            search = self.idx.search([3014687], [3014690])
            self.assertEqual(len(list(search)), 2)
            search = self.idx.search([3014687], [3014691])
            self.assertEqual(len(list(search)), 2)

        def test_correct_block_length(self):
            """Following issues 504 and 1086.

            https://github.com/biopython/biopython/pull/504
            https://github.com/biopython/biopython/pull/1086#issuecomment-285080702

            We get the alignement corresponding to the following whole MAF block
            and check that the lengths of its sequences are correct:

            a score=40840.000000
            s mm9.chr10                         3014689 53 + 129993255 GGGAGCATAAAACTCTAAATCTGCTAAATGTCTTGTCCCT-TTGGAAAGAGTTG
            s hg18.chr6                        15870832 53 - 170899992 GGGATCATAAACCATTTAATCTGTGAAATATCTAATCTTT-TGGGAAATAGTGG
            i hg18.chr6                        C 0 I 401
            s panTro2.chr6                     16389401 53 - 173908612 GGGATCATAAACCATTTAATCTGTGAAATATCTAATCTTT-TGGGAAATAGTGG
            q panTro2.chr6                                             9999999999999999999999999999999999999999-9999999999999
            i panTro2.chr6                     C 0 I 400
            s calJac1.Contig6394                   6228 53 +    133105 GGGATCATAAGCCATTTAATCTGTGAAATGTGAAATCTTT-TGGGAAACAGTGG
            i calJac1.Contig6394               C 0 I 2
            s otoGar1.scaffold_334.1-359464      184148 52 -    359464 GGAAGCATAAACT-TTTAATCTATGAAATATCAAATCACT-TGGGCAATAGCTG
            q otoGar1.scaffold_334.1-359464                            7455455669566-99665699769895555689997599-9984787795599
            i otoGar1.scaffold_334.1-359464    I 2931 I 2
            s loxAfr1.scaffold_75566               1201 54 -     10574 GGGAGTATAAACCATTTAGTCTGCGAAATGCCAAATCTTCAGGGGAAAAAGCTG
            q loxAfr1.scaffold_75566                                   899989799999979999999999999999797999999999999999999999
            i loxAfr1.scaffold_75566           C 0 I 2
            e tupBel1.scaffold_114895.1-498454   167376 4145 -    498454 I
            e echTel1.scaffold_288249             87661 7564 +    100002 I
            e ponAbe2.chr6                     16161448 8044 - 174210431 I
            """
            ali = self.idx.get_spliced([3014689], [3014689 + 53])
            seq_dict = dict([(seqrec.id, seqrec.seq) for seqrec in ali])
            correct_lengths = {
                "mm9.chr10": 53,
                "hg18.chr6": 53,
                "panTro2.chr6": 53,
                "calJac1.Contig6394": 53,
                "otoGar1.scaffold_334.1-359464": 52,
                "loxAfr1.scaffold_75566": 54}
            for seq_id, length in correct_lengths.items():
                self.assertEqual(len(seq_dict[seq_id].ungap('-')), length)

        def test_correct_spliced_sequences_1(self):
            """Checking that spliced sequences are correct.

            We get the alignement corresponding to the following whole MAF block
            and check that the sequences are correct:

            a score=40840.000000
            s mm9.chr10                         3014689 53 + 129993255 GGGAGCATAAAACTCTAAATCTGCTAAATGTCTTGTCCCT-TTGGAAAGAGTTG
            s hg18.chr6                        15870832 53 - 170899992 GGGATCATAAACCATTTAATCTGTGAAATATCTAATCTTT-TGGGAAATAGTGG
            i hg18.chr6                        C 0 I 401
            s panTro2.chr6                     16389401 53 - 173908612 GGGATCATAAACCATTTAATCTGTGAAATATCTAATCTTT-TGGGAAATAGTGG
            q panTro2.chr6                                             9999999999999999999999999999999999999999-9999999999999
            i panTro2.chr6                     C 0 I 400
            s calJac1.Contig6394                   6228 53 +    133105 GGGATCATAAGCCATTTAATCTGTGAAATGTGAAATCTTT-TGGGAAACAGTGG
            i calJac1.Contig6394               C 0 I 2
            s otoGar1.scaffold_334.1-359464      184148 52 -    359464 GGAAGCATAAACT-TTTAATCTATGAAATATCAAATCACT-TGGGCAATAGCTG
            q otoGar1.scaffold_334.1-359464                            7455455669566-99665699769895555689997599-9984787795599
            i otoGar1.scaffold_334.1-359464    I 2931 I 2
            s loxAfr1.scaffold_75566               1201 54 -     10574 GGGAGTATAAACCATTTAGTCTGCGAAATGCCAAATCTTCAGGGGAAAAAGCTG
            q loxAfr1.scaffold_75566                                   899989799999979999999999999999797999999999999999999999
            i loxAfr1.scaffold_75566           C 0 I 2