Python pysam.asGTF函数代码示例

 def testDisjointIterators(self):
     # two iterators working on the same file
     tabix = pysam.TabixFile(self.filename)
     a = tabix.fetch(parser=pysam.asGTF(), multiple_iterators=True).next()
     b = tabix.fetch(parser=pysam.asGTF(), multiple_iterators=True).next()
     # both iterators are at top of file
     self.assertEqual(str(a), str(b))

    def testJoinedIterators(self):

        # two iterators working on the same file
        tabix = pysam.TabixFile(self.filename)
        a = tabix.fetch(parser=pysam.asGTF()).next()
        b = tabix.fetch(parser=pysam.asGTF()).next()
        # the first two lines differ only by the feature field
        self.assertEqual(a.feature, "UTR")
        self.assertEqual(b.feature, "exon")
        self.assertEqual(re.sub("UTR", "", str(a)), re.sub("exon", "", str(b)))

def _TestMultipleIteratorsHelper(filename, multiple_iterators):
    """open file within scope, return iterator."""

    tabix = pysam.TabixFile(filename)
    iterator = tabix.fetch(parser=pysam.asGTF(), multiple_iterators=multiple_iterators)
    tabix.close()
    return iterator

def read_gtf(lines, scaffolds, contig_prefix):
    table = {}  # gene_id -> transcript_id -> exon_number -> feature -> [items]
    for gtf in text.parse_lines(lines, pysam.asGTF()):
        if not filter_gtf_record(gtf):
            update_gtf_table(table, gtf, scaffolds, contig_prefix)

    return table

 def subset_of_feature_in_region(self, contig=None, start=None, end=None, types=None):
     '''
     Example:
      # return dict of rna and tss id
      for rna_tss in tabix.subset_of_feature_in_region(contig="chr2L", end=9839,
                                                 types=["transcript_id", 'tss_id']):
          print rna_tss
     '''
     
     for gtf in pysam.Tabixfile.fetch(self.tabixfile, contig, start, end,
                                      parser=pysam.asGTF()):
         if isinstance(types, str):
             try:
                 yield gtf.asDict()[types]
             except KeyError:
                 print 'key \'{0}\' is not found in {1}'.format(t, self.ingtf)
                 
         elif isinstance(types, list):
             tmp = dict()
             for t in types:
                 try:
                     tmp.update({t: gtf.asDict()[t]})
                 except KeyError:
                     print 'key \'{0}\' is not found in {1}'.format(t, self.ingtf)
             yield tmp

 def gene_in_region(self, contig=None, start=None, end=None):
     for gtf in pysam.Tabixfile.fetch(self.tabixfile, contig, start, end,
                                      parser=pysam.asGTF()):
         try:
             yield gtf.asDict()['gene_name']
         except KeyError:
             print 'key \'{0}\' is not found in {1}'.format(t, self.ingtf)

    def testCopy(self):
        a = self.tabix.fetch(parser=pysam.asTuple()).next()
        b = copy.copy(a)
        self.assertEqual(a, b)

        a = self.tabix.fetch(parser=pysam.asGTF()).next()
        b = copy.copy(a)
        self.assertEqual(a, b)

    def __init__(self, file_path):
        file_path = str(file_path)
        if not file_path.endswith('.gz'):
            if os.path.exists(file_path + '.gz'):
                file_path += '.gz'
            else:
                file_path = self.compress(file_path)

        super().__init__(file_path, parser=pysam.asGTF())

    def _open_file(self):  # type: (...) -> pysam.TabixFile
        # Open gtf file.
        gtf_file = pysam.TabixFile(
            native_str(self._gtf_path), parser=pysam.asGTF())

        # Yield file object and ensure it is closed.
        try:
            yield gtf_file
        finally:
            gtf_file.close()

def load_kg_gtf(gtf_file_name):
    f = pysam.TabixFile(gtf_file_name)
    gtf = f.fetch(parser=pysam.asGTF())
    feats = []
    for row in gtf:
        attr = parse_gtf_attr(row.attributes)
        currfeat = GFFFeature(row.seqname, row.source, row.feature,\
                                int(row.start), int(row.end), score, \
                                strand, frame, attr)
        feats.append(currfeat)
    return feats

    def testRead(self):

        for x, r in enumerate(self.tabix.fetch(parser=pysam.asGTF())):
            c = self.compare[x]
            self.assertEqual(len(c), len(r))
            self.assertEqual(list(c), list(r))
            self.assertEqual(c, str(r).split("\t"))
            self.assertTrue(r.gene_id.startswith("ENSG"))
            if r.feature != 'gene':
                self.assertTrue(r.transcript_id.startswith("ENST"))
            self.assertEqual(c[0], r.contig)

    def testSetting(self):

        for r in self.tabix.fetch(parser=pysam.asGTF()):
            r.contig = r.contig + "_test"          
            r.source = r.source + "_test"
            r.feature = r.feature + "_test"
            r.start += 10
            r.end += 10
            r.score = 20
            r.strand = "+"
            r.frame = 0
            r.attributes = 'gene_id "0001";'

 def fetch_gtf(self, contig=None, start=None, end=None):
     '''
     Returns:
      pysam.tabix object
     
     Example:
      # return dict of each GTF line
      for _ in tabix.fetch_gtf(contig="chr2L", end=9839):
          print _.asDict()
     '''
     
     for gtf in pysam.Tabixfile.fetch(self.tabixfile, contig, start, end,
                                      parser=pysam.asGTF()):
         yield gtf

def overlap_annotation(junc, anno):
    ts_set = anno.fetch(junc.chrom, junc.start, junc.end, parser=pysam.asGTF())
    exon_set = filter(lambda x: x.feature=='exon', [i for i in ts_set])
    for idx,exon in enumerate(exon_set):
        if idx == len(exon_set) - 1:
            continue
        

    for(std::size_t i = 0; i < exons.size(); i++) {
        if(exons[i].start > junction.end) {
            //No need to look any further
            //the rest of the exons are outside the junction
            break;
        }
        //known junction
        if(exons[i].end == junction.start &&
                exons[i + 1].start == junction.end) {
            junction.known_acceptor = true;
            junction.known_donor = true;
            junction.known_junction = true;
            known_junction = true;
        }
        else {
            if(!junction_start) {
                if(exons[i].end >= junction.start) {
                    junction_start = true;
                }
            }
            if(junction_start) {
                if(exons[i].start > junction.start &&
                        exons[i].end < junction.end) {
                    junction.exons_skipped.insert(exons[i].name);
                }
                if(exons[i].start > junction.start) {
                    junction.donors_skipped.insert(exons[i].start);
                }
                if(exons[i].end < junction.end) {
                    junction.acceptors_skipped.insert(exons[i].end);
                }
                if(exons[i].end == junction.start) {
                    junction.known_donor = true;
                }
                //TODO - check for last exon
                if(exons[i].start == junction.end) {
                    junction.known_acceptor = true;
                }
            }
        }
    }

    def from_gtf(
            cls,
            gtf_path,  # type: pathlib.Path
            chromosomes=None,  # type: List[str]
            record_filter=None  # type: Callable[[Any], bool]
    ):  # type: (...) -> TranscriptReference
        """Builds an Reference instance from the given GTF file."""

        # Open gtf file.
        gtf = pysam.TabixFile(native_str(gtf_path), parser=pysam.asGTF())

        if chromosomes is None:
            chromosomes = gtf.contigs

        # Build the trees.
        transcript_trees = {}
        exon_trees = {}

        for chrom in chromosomes:
            # Collect exons and transcripts.
            transcripts = []
            exons = []

            records = gtf.fetch(reference=chrom)

            if record_filter is not None:
                records = (rec for rec in records if record_filter(rec))

            for record in records:
                if record.feature == 'transcript':
                    transcripts.append(cls._record_to_transcript(record))
                elif record.feature == 'exon':
                    exons.append(cls._record_to_exon(record))

            # Build transcript lookup tree.
            transcript_trees[chrom] = IntervalTree.from_tuples(
                (tr.start, tr.end, tr) for tr in transcripts)

            # Build exon lookup tree.
            keyfunc = lambda rec: rec.transcript_id

            exons = sorted(exons, key=keyfunc)
            grouped = itertools.groupby(exons, key=keyfunc)

            for tr_id, grp in grouped:
                exon_trees[tr_id] = IntervalTree.from_tuples(
                    (exon.start, exon.end, exon) for exon in grp)

        return cls(transcript_trees, exon_trees)