专栏首页生信菜鸟团python:手动比对序列并绘制测序饱和度图片

python:手动比对序列并绘制测序饱和度图片

最近因为工作需要,有一组RNA探针测序数据要求检查其测序饱和度的情况,来评估测序的冗余度。

测序饱和度的评估参考RNA-seq的定义,并非10X定义的根据UMI计算的测序饱和度。

3)饱和度评估 得到sam文件后其实我们可以计算出每一条reads落在哪里,基因间区或者基因区(基因区的话在哪个基因上也可以列出来)。假如我们得到3,200,000条uniq mapping的reads。我们可以采用梯度为10,000,分别为100,000, 200,000, 300,000, 。。。3,100,000, 3,200,000,依次随即抽样,看抽出来的这些reads分别检测到多少基因。然后把reads数和检测的基因数画一个曲线,看看这条曲线是否饱和。若,当reads数到一定程度检测到的基因几乎不增加或者很少增加,则测序包和,否者就是测序量不够没有达到饱和。

由于测序数据是探针数据,并且数量也不是太多,考虑使用python的正则进行序列匹配,实际结果看其比对效率还是挺低的。具体如下:

只以一个测序文件为例,gzip包用于对fastq.gz格式的测序文件进行读取处理。

import gzip
import os
import re
import random

# fastq.gz文件路径,只以一个测序文件为例
os.chdir("F:\\python\\测序饱和度")
fastq_file = [i for i in os.listdir() if i.endswith('gz')]
fastq_file = fastq_file[0]

处理模板序列

模板序列都是12bp的短序列,共有8684个唯一探针序列。每两行是一个探针信息,第一行是以">"开头的探针名称,第二行是具体序列。探针的两行信息以"|分隔合并为一个字符串。为了提高正则匹配的效率,将所有的模板探针序列以逗号分隔并成一个字符串。

# 示例探针
# >RNA0048103
# AAGTGCTCCTGA

ref_seq = []
r = ""
with open('rna.fa') as f:
    for line in f:
        if line.startswith('>'):
            r = line.strip()
        else:
            ref_seq.append("{}|{}".format(r, line.strip()))
            r = ""

ref_seq_str = ",".join(ref_seq)

# 模板序列字符串
ref_seq_str[:50]
#'>RNA0048103|AAGTGCTCCTGA,>RNA0048199|TGCGGAGTTAGA,'

提取测序序列

使用循环处理fastq.gz文件,提取序列至一个列表中做后续处理。

使用n来控制提取序列,遇到@开头的行,则将n标记为1,下一次循环时则提取整行数据,将其置于预先定义的列表中。

# 示例
# @TPNB500160AR:324:H5VCKBGXJ:1:11101:10065:1093 2:N:0:CGTGGCTT+TTAATGAT
# GTGTCACTGTATCGGGCAATCACAAAA
# +
# AAAAAAEEEAEEEEEEEEEEEEEEEE/

fastq_seq = []
n = 0
if(os.path.exists(fastq_file)):
    with gzip.open(fastq_file) as gf:
        for line in gf:
            if line.startswith(b'@'): 
                n = 1
                continue
            if n == 1:
                fastq_seq.append(line.strip().decode()[14:26])
                n = 0
                continue
else:
    print("File {} not exitst!".format(fastq_file))

正则进行序列比对

使用正则进行序列匹配,如果匹配,则返回探针序号,如果没有匹配,则返回字符串“None”。

共有86完条read,比对共运行接近8min,效率比较低,使用常规字符串操作进行序列匹配还是只适用于数据量比较少的情况。

def align(fastq_seq, ref_seq):
    s = []
    for seq in fastq_seq:
        pat = "(RNA[0-9]+)\|{}".format(seq)
        rr = re.search(pat, ref_seq)
        if rr:
            s.append(rr.group(1))
        else:
            s.append("None")
    return s

res_align = align(fastq_seq, ref_seq_str)

res_align[:3]
# ['RNA42588',
# 'RNA46325',
# 'None']

随机抽样获取饱和度数据

饱和度数据其实就是重抽样数据,筛选到比对结果后,去重并计数。

dat = []
# 抽样数量1000 - 10000,1000为步长
for i in range(100000//1000 + 1):
    res = random.sample(res_align, i * 1000)
    l = list(set(res)) # 去重
    l = [i for i in l if i != "None"] # 去除比对失败的探针
    dat.append([i * 1000, len(l)]) # 

# 抽样数量10,0000 - 总reads数,100000为步长
fastq_seq_len = len(fastq_seq) // 100000
for i in range(fastq_seq_len + 1):
    res = random.sample(res_align, i * 100000)
    l = list(set(res))
    l = [i for i in l if i != "None"]
    dat.append([i * 100000, len(l)])

# 查看重抽样的结果
dat[:3]
# [[0, 0], [1000, 784], [2000, 1399]]

绘制散点图

import pandas as pd
import matplotlib.pyplot as plt

pd_dat = pd.DataFrame(dat, columns=['Reads number', 'Detected probe'])
pd_dat["Reads number"] = pd_dat["Reads number"]/10000

plt.scatter(x=pd_dat["Reads number"], y=pd_dat["Detected probe"])
plt.xlabel("Reads number (X10000)")
plt.ylabel("Detected probe")

本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:冰糖

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-07-27

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用Python+opencv进行图像处理(一) | 视觉入门

    计算机视觉是人工智能最热门的应用领域之一。人工智能技术推动了汽车自动驾驶、机器人以及各种照片处理类软件的巨大发展。目标检测技术也在稳步推进。生成对抗网络(GAN...

    磐创AI
  • 解密ATAC中的测序饱和度分析

    测序数据量对于NGS数据分析是非常重要的,测序数据量过低,不能有效覆盖基因组完整信息,测序数据量过高,则会造成冗余,不够经济。为了验证当前测序量能否满足需求,或...

    生信修炼手册
  • 这10个开源人工智能项目,你必须了解!

    来自:开源中国 链接:https://my.oschina.net/editorial-story/blog/1592254 推荐 10 个饱受好评且功能独特的...

    企鹅号小编
  • 使用颜色空间进行图像分割

    原文地址:https://realpython.com/python-opencv-color-spaces/

    努力努力再努力F
  • 开发者不可错过,10个简单的开源项目——人工智能篇

    AI UNION 人工智能产业技术创新战略联盟 这里是人工智能联盟,汇聚了最新的AI新闻资讯,还有最前沿的国内外AI开源技术,最具价值的AI创新企业,最具权威的...

    企鹅号小编
  • 开发者不可错过的开源项目 —— 人工智能篇

    推荐 10 个饱受好评且功能独特的开源人工智能项目 关于人工智能的项目,相信大家都看过或者用过不少了,但它们的大多数看上去都十分“高大上”,让人感觉要掌握他们犹...

    CSDN技术头条
  • cell ranger分析结果详细解读

    输出文件非常的多,为了方便查看结果,提供了一个所有结果汇总的html页面,即web_summary.html。该网页的结果分成了summary和analysis...

    生信修炼手册
  • .NET3.5 GDI+ 图形操作1

          前言: 本文章抄袭自本人刚刚买的《ASP.NET 3.5从入门到精通》这本书,此书介绍在 http://www.china-pub.com/4499...

    跟着阿笨一起玩NET
  • 转录组分析的正确知识都了解了吗?

    转录组分析是目前应用最广的高通量测序分析技术之一。常见设计是不同样品之间比较,寻找差异基因、标志基因、协同变化基因、差异剪接和新转录本,并进行结果可视化、功能注...

    生信宝典
  • 数据可视化的基本流程总结

    我们要的不是数据,而是数据告诉我们的事实。大多数人面临这样一个挑战:我们认识到数据可视化的必要性,但缺乏数据可视化方面的专业技能。部分原因可以归结于,数据可视化...

    1480
  • 简单实现Android刮刮卡效果

    本文实例为大家分享了Android仿刮刮卡效果展示的具体代码,供大家参考,具体内容如下

    砸漏
  • 数据可视化(11)-Seaborn系列 | 小提琴图violinplot()

    小提琴形图(violin plot)的作用与盒形图(box plot)和whidker plot的作用类似,它显示了一个或多个分类变量的几个级别的定量数据的分布...

    数据分析可视化
  • 10行代码实现python人脸识别

    人脸识别,是基于人的脸部特征信息进行身份识别的一种生物识别技术。用摄像机或摄像头采集含有人脸的图像或视频流,并自动在图像中检测和跟踪人脸,进而对检测到的人脸进行...

    星星在线
  • 数据可视化(5)-Seaborn系列 | 柱状图countplot()

    案例代码已上传:Github https://github.com/Vambooo/SeabornCN

    数据分析可视化
  • 转录组分析的正确姿势(第三版)

    转录组分析是目前应用最广的高通量测序分析技术之一。常见设计是不同样品之间比较,寻找差异基因、标志基因、协同变化基因、差异剪接和新转录本,并进行结果可视化、功能注...

    生信宝典
  • 转录组分析的正确姿势

    生信宝典
  • 5秒钟内将手绘网站线框图转换为可用的 HTML网站

    你可以在 GitHub 上找到这个项目的代码:https://github.com/ashnkumar/sketch-code

    IT派
  • 资源 | 深度学习自动前端开发:从草图到HTML只需5秒(附代码)

    选自InsightDataScience 作者:Ashwin Kumar 机器之心编译 参与:乾树、李泽南 在人们的不断探索下,「使用人工智能自动生成网页」的方...

    机器之心
  • 我们整理了20个Python项目,送给正在求职的你

    职场中一贯有“金三银四”、“金九银十”的说法。如果你是一名正在求职或准备跳槽的程序员,不妨趁着这两个月时间好好准备一下。

    用户7886150

扫码关注云+社区

领取腾讯云代金券