美文网首页
HTSeq中GenomicInterval和GenomicPos

HTSeq中GenomicInterval和GenomicPos

作者: scdzzdw | 来源:发表于2019-10-28 15:20 被阅读0次

    举个栗子,我有三个位点,每个位点赋值为2,但是位点有重叠,如下

    import HTSeq
    
    ga = HTSeq.GenomicArray("auto", stranded=False)
    
    iv = HTSeq.GenomicInterval('chr1', 50, 51, "+")
    ga[iv] += 2
    iv = HTSeq.GenomicInterval('chr1', 50, 51, "+")
    ga[iv] += 2
    iv = HTSeq.GenomicInterval('chr1', 50, 52, "+")
    ga[iv] += 2
    iv = HTSeq.GenomicInterval('chr1', 49, 53, "+")
    ga[iv] += 2
    
    for iv, value in ga.steps():
        if value:
            print(iv.chrom, iv.start, iv.end, value)
    

    得出地结果为

    chr1 49 50 2.0
    chr1 50 51 8.0
    chr1 51 52 4.0
    chr1 52 53 2.0
    

    但是我想要的输出结果如下

    chr1 49 53 2.0
    chr1 50 51 4.0
    chr1 50 52 2.0
    

    使用GenomicPosition进行位置储存,代码如下

    ga = HTSeq.GenomicArray("auto", stranded=False)
    iv = HTSeq.GenomicPosition('chr1', 50, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 50, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 49, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 53, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 48, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 52, "+")
    ga[iv] += 1
    
    for iv, value in ga.steps():
        if value:
            print(iv, iv.start, iv.end, value)
    

    输出结果为

    chr1:[48,50)/. 48 50 1.0
    chr1:[50,51)/. 50 51 2.0
    chr1:[52,54)/. 52 54 1.0
    

    可以看到,position 49我赋予了值1,但是chr1:[48,50)/. 48 50 1.0把49覆盖了,那么

    如何修改?

    ============================================
    尝试了两种方法,最后只能通过每个位点进行记录保存

    ga = HTSeq.GenomicArray("auto", stranded=False)
    iv = HTSeq.GenomicPosition('chr1', 50, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 50, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 49, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 53, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 48, "+")
    ga[iv] += 1
    iv = HTSeq.GenomicPosition('chr1', 52, "+")
    ga[iv] += 1
    
    for iv, value in ga.steps():
        if value:
            for i in range(iv.start, iv.end):
                print(i, value)
    

    相关文章

      网友评论

          本文标题:HTSeq中GenomicInterval和GenomicPos

          本文链接:https://www.haomeiwen.com/subject/wiuvvctx.html