使用 python 进行 Ribo-seq 质控分析 二

测试开头


接纳自己,放过自己。

1引言

这部分是绘制 metagene plot,参考前面那篇 cell 文章的方法,并对代码进行改写和优化,使之适用于自己的数据。

2from to start codon

定义函数:

def metagenePfofileToST(inputFile,outputFile,cdslength=400,expression=50):
    #########################################################
    # prepare gene expression and length
    with open(inputFile,'r') as input:
        gene_infoDict = {}
        for line in input:
            fileds = line.split()
            gene_name = fileds[0].split('_')[0]
            pos = int(fileds[1])
            start = int(fileds[0].split('_')[3])
            end = int(fileds[0].split('_')[4])
            cdsLength = end - start + 1
            # geneLength = int(fileds[0].split('|')[4])
            density = float(fileds[2])
            # filter CDS > 400 nt gene
            if cdsLength > cdslength:
                # key
                key = ':'.join([gene_name,str(cdsLength)])
                gene_infoDict.setdefault(key,0)
                if start <= pos <= end:
                    gene_infoDict[key] += density
                else:
                    pass
            else:
                pass

    # filter CDS expression > 50
    filtedGeneDict = {}
    for key,val in gene_infoDict.items():
        if val > expression:
            meanNorm = val/int(key.split(':')[1])
            filtedGeneDict[key] = [val,meanNorm]
        else:
            pass
    #########################################################
    # Meta-gene analysis from start codon
    mylist = range(-50, 1501)
    rangeDict = dict([i, 0] for i in mylist)
    countDict = dict([i, 0] for i in mylist)

    # open file
    with open(inputFile,'r') as input:
        gene_infoDict = {}
        for line in input:
            fileds = line.split()
            gene_name = fileds[0].split('_')[0]
            pos = int(fileds[1])
            start = int(fileds[0].split('_')[3])
            end = int(fileds[0].split('_')[4])
            cdsLength = end - start + 1
            # geneLength = int(fileds[0].split('|')[4])
            density = float(fileds[2])
            id = ':'.join([gene_name,str(cdsLength)])
            if id in filtedGeneDict:
                reldist = pos - start
                if -50 <= reldist <= 1500:
                    # divide reads at one position by average number of reads per base for this gene
                    reads = density / filtedGeneDict[id][1]
                    rangeDict[reldist] += reads
                    # how often was this position counted in the calculation
                    countDict[reldist] += 1
                else:
                    pass
            else:
                pass

    #########################################################
    # output data
    # sort dict
    tupledlist1 = list(rangeDict.items())
    tupledlist1.sort()
    tupledlist2 = list(countDict.items())
    tupledlist2.sort()

    fullDict = {}
    zippedlist = zip(tupledlist1, tupledlist2)
    for elem in zippedlist:
        col0 = elem[0][0]       # list0 col0 = position (K)
        col1 = elem[0][1]       # list0 col1 = norm read number
        col2 = elem[1][1]       # list1 col1 = how often was position counted

        #normalization2 by frequnecy
        if col2 == 0:
            fullDict[col0] = 0
        else:
            fullDict[col0] = col1 / col2

    # Finish output
    tupledlist = list(fullDict.items())
    tupledlist.sort()

    # output
    outFileP = open(outputFile, 'w')

    for elem in tupledlist:
        outFileP.write('t'.join([str(elem[0]),str(elem[1])]) + 'n')
    outFileP.close()

分析:

metagenePfofileToST('ribo-EB.density.txt','ribo-EB.metegene2StartCodon.txt',cdslength=600,expression=100)

metagenePfofileToST('ribo-ESC.density.txt','ribo-ESC.metegene2StartCodon.txt',cdslength=600,expression=100)

3from to start codon

定义函数:

def metagenePfofileToSP(inputFile,outputFile,cdslength=400,expression=50):
    #########################################################
    # prepare gene expression and length
    with open(inputFile,'r') as input:
        gene_infoDict = {}
        for line in input:
            fileds = line.split()
            gene_name = fileds[0].split('_')[0]
            pos = int(fileds[1])
            start = int(fileds[0].split('_')[3])
            end = int(fileds[0].split('_')[4])
            cdsLength = end - start + 1
            # geneLength = int(fileds[0].split('|')[4])
            density = float(fileds[2])
            # filter CDS > 400 nt gene
            if cdsLength > cdslength:
                # key
                key = ':'.join([gene_name,str(cdsLength)])
                gene_infoDict.setdefault(key,0)
                if start <= pos <= end:
                    gene_infoDict[key] += density
                else:
                    pass
            else:
                pass

    # filter CDS expression > 50
    filtedGeneDict = {}
    for key,val in gene_infoDict.items():
        if val > expression:
            meanNorm = val/int(key.split(':')[1])
            filtedGeneDict[key] = [val,meanNorm]
        else:
            pass
    #########################################################
    # Meta-gene analysis from stop codon
    mylist = range(-1500, 51)
    rangeDict = dict([i, 0] for i in mylist)
    countDict = dict([i, 0] for i in mylist)

    # open file
    with open(inputFile,'r') as input:
        gene_infoDict = {}
        for line in input:
            fileds = line.split()
            gene_name = fileds[0].split('_')[0]
            pos = int(fileds[1])
            start = int(fileds[0].split('_')[3])
            end = int(fileds[0].split('_')[4])
            cdsLength = end - start + 1
            # geneLength = int(fileds[0].split('|')[4])
            density = float(fileds[2])
            id = ':'.join([gene_name,str(cdsLength)])
            if id in filtedGeneDict:
                # reldist = pos - start
                reldisp = pos - (end + 1)
                if -1500 <= reldisp <= 50:
                    # divide reads at one position by average number of reads per base for this gene
                    reads = density / filtedGeneDict[id][1]
                    rangeDict[reldisp] += reads
                    # how often was this position counted in the calculation
                    countDict[reldisp] += 1
                else:
                    pass
            else:
                pass

    #########################################################
    # output data
    # sort dict
    tupledlist1 = list(rangeDict.items())
    tupledlist1.sort()
    tupledlist2 = list(countDict.items())
    tupledlist2.sort()

    fullDict = {}
    zippedlist = zip(tupledlist1, tupledlist2)
    for elem in zippedlist:
        col0 = elem[0][0]       # list0 col0 = position (K)
        col1 = elem[0][1]       # list0 col1 = norm read number
        col2 = elem[1][1]       # list1 col1 = how often was position counted

        #normalization2 by frequnecy
        if col2 == 0:
            fullDict[col0] = 0
        else:
            fullDict[col0] = col1 / col2

    # Finish output
    tupledlist = list(fullDict.items())
    tupledlist.sort()

    # output
    outFileP = open(outputFile, 'w')

    for elem in tupledlist:
        outFileP.write('t'.join([str(elem[0]),str(elem[1])]) + 'n')
    outFileP.close()

分析:

metagenePfofileToSP('ribo-EB.density.txt','ribo-EB.metegene2StopCodon.txt',cdslength=600,expression=100)

metagenePfofileToSP('ribo-ESC.density.txt','ribo-ESC.metegene2StopCodon.txt',cdslength=600,expression=100)

4绘图

批量导入数据:

library(ggplot2)
library(tidyverse)
library(ggsci)
library(Rmisc)
library(data.table)

###################################################################### start codon
name <- list.files(path = './',pattern = 'metegene2StartCodon.txt')

map_df(1:length(name), function(x){
  tmp = read.table(paste('./',name[x],sep = ''))
  colnames(tmp) <- c('pos','density')
  tmp$sample <- sapply(strsplit(name[x],split = '\.'),'[',1)
  return(tmp)
}) %>% data.table() -> df_st

# check
head(df_st,3)
#    pos  density  sample
# 1: -50 1.551956 ribo-EB
# 2: -49 1.589671 ribo-EB
# 3: -48 1.581825 ribo-EB

绘图:

# single group
ggplot(df_st,aes(x = pos,y = density)) +
  geom_line(aes(color = sample)) +
  geom_hline(yintercept = 1,lty = 'dashed',size = 1,color = 'grey50') +
  theme_classic(base_size = 16) +
  # scale_color_d3(name = '') +
  scale_color_futurama(name = '') +
  theme(legend.position = c(0.88,0.88),
        legend.background = element_blank(),
        aspect.ratio = 0.6,
        strip.background = element_rect(color = NA,fill = 'grey')) +
  # facet_grid(exp~type) +
  ylab('Mean read density [AU]') +
  xlab('Distance from start codon (nt)')
使用 python 进行 Ribo-seq 质控分析 二

对于 from to stop codon 也是类似:

###################################################################### stop codon
name <- list.files(path = './',pattern = 'metegene2StopCodon.txt')

map_df(1:length(name), function(x){
  tmp = read.table(paste('./',name[x],sep = ''))
  colnames(tmp) <- c('pos','density')
  tmp$sample <- sapply(strsplit(name[x],split = '\.'),'[',1)
  return(tmp)
}) %>% data.table() -> df_sp

# check
head(df_sp,3)
#      pos  density  sample
# 1: -1500 1.512251 ribo-EB
# 2: -1499 1.608943 ribo-EB
# 3: -1498 1.499345 ribo-EB

# single group
ggplot(df_sp,aes(x = pos,y = density)) +
  geom_line(aes(color = sample)) +
  geom_hline(yintercept = 1,lty = 'dashed',size = 1,color = 'grey50') +
  theme_classic(base_size = 16) +
  # scale_color_d3(name = '') +
  scale_color_futurama(name = '') +
  theme(legend.position = c(0.88,0.88),
        legend.background = element_blank(),
        aspect.ratio = 0.6,
        strip.background = element_rect(color = NA,fill = 'grey')) +
  # facet_grid(exp~type) +
  ylab('Mean read density [AU]') +
  xlab('Distance from stop codon (nt)')
使用 python 进行 Ribo-seq 质控分析 二

不知道为啥有两个特别高的峰,感觉还挺奇怪。

5结尾

质控部分差不多就到这里结束了。



使用 python 进行 Ribo-seq 质控分析 二


欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。

群二维码：

使用 python 进行 Ribo-seq 质控分析 二


老俊俊微信：


使用 python 进行 Ribo-seq 质控分析 二

知识星球：


使用 python 进行 Ribo-seq 质控分析 二


所以今天你学习了吗？

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！



 往期回顾 





◀deeptools 结果重新绘图

◀使用 python 进行 Ribo-seq 质控分析 一

◀跟着 Cell 学 Ribo-seq 分析 四 (终)

◀跟着 Cell 学 Ribo–seq 分析 三 (Metagene Plot)

◀跟着 Cell 学 Ribo–seq 分析 二

◀RiboPlotR 优雅的可视化你的 Ribo–seq 数据

◀跟着 Cell 学 Ribo–seq 分析 一

◀RNAmod: m6A peak 注释神器

◀提取酵母两端扩展 50nt 的 CDS 序列

◀R 里使用 python 加速 R 代码运行

◀...

测试结尾

请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！

作者头像
weinfoadmin

你可能也喜欢

8-1651542331
跟着Nature学绘图(2) 箱线图-累积分布曲线图
2 5月, 2022
9-1651542322
Julia 笔记之字符串
2 5月, 2022
0-1651542343
Julia 笔记之数学运算和初等函数
1 5月, 2022