生信分析进阶2 - 利用GC含量的Loess回归矫正reads数量

在NGS数据比对后,需要矫正GC偏好引起的reads数量误差可用loess回归算法,使用R语言对封装的loess算法实现。

在NIPT中,GC矫正对检测结果准确性非常重要,具体研究参考以下文章。

Noninvasive Prenatal Diagnosis of Fetal Trisomy 18 and Trisomy 13 by Maternal Plasma DNA Sequencing
链接地址:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3130771/
在这里插入图片描述窗口划分可参考文章:

生信软件8 - bedtools进行窗口划分、窗口GC含量、窗口测序深度和窗口SNP统计

获取参考基因组大小

以hg19参考基因组为例。

# 安装python库
pip install pyfaidx# 保留chr1-chr22 chrX chrY
faidx reference/hg19.fasta -i chromsizes|grep -E -v '_|chrM' > hg19.genome.size

hg19.genome.size

划分基因组窗口

以1000kb划分为例。

bedtools makewindows -g hg19.genome.size -w 1000000 > hg19.1000kb.bed

hg19.1000kb.bed

划分窗口

bedtools nuc -fi /reference/hg19.fasta -bed hg19.1000kb.bed|cut -f 1-3,5 > hg19.1000kb.gc.bed

hg19.1000kb.gc.bed

统计窗口reads和GC含量

bedtools coverage -a hg19.1000kb.bed -b sample.sorted.bam > sample.count

sample.count

整理数据

paste <(grep -v '#' hg19.1000kb.gc.bed) <(cut -f4 sample.count)|sed '1i chr\tstart\tend\tGC\treads' > sample.gc.count

sample.gc.count

利用GC含量的Loess回归矫正reads数量

R语言实现。

# loess_gc_correct.R
# Useage: Rscript loess_gc_correct.R /path/sample.gc.count /path/sample.gc.corrected.countargs=commandArgs(T)
# 输入文件路径
gc.reads.file <- args[1]
# 输出文件路径
gc.reads.corrected.file <- args[2]# 读取输入文件
raw.data <- read.table(gc.reads.file, sep='\t', head=TRUE)# loess regression 进行GC矫正reads数量
gc.count.loess <- loess(reads~GC,data = raw.data,control = loess.control(surface = "direct"), degree=2) prediction <- predict(gc.count.loess, raw.data$GC)raw.data$corrected_reads <- as.integer(prediction)# 保存
write.table(raw.data, file = gc.reads.corrected.file,sep = '\t', quote = FALSE)

矫正后结果

矫正后文件

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/323272.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

如何把多个文件(夹)平均复制到多个文件夹中

首先&#xff0c;需要用到的这个工具&#xff1a; 度娘网盘 提取码&#xff1a;qwu2 蓝奏云 提取码&#xff1a;2r1z 假定的情况是&#xff0c;共有20个兔兔的图片&#xff0c;想要平均的复制4个文件夹里&#xff0c;那么每个文件夹里面就有5个图片 &#xff08;如果是5个&a…

nginx自动部署-跨操作系统

项目里面有一个需求&#xff0c;就是需要用让nginx进程提供给系统管理一个start,stop和getPid方法&#xff0c;这样系统管理可以自动拉起来nginx&#xff0c;达到自动部署的目的。离线部署同样适用 这样一来&#xff0c;我就需要提供windows版本linux不同版本的nginx源码包&am…

Spring JdbcTemplate实现自定义动态sql拼接功能

需求描述&#xff1a; sql 需要能满足支持动态拼接&#xff0c;包含 查询字段、查询表、关联表、查询条件、关联表的查询条件、排序、分组、去重等 实现步骤&#xff1a; 1&#xff0c;创建表及导入测试数据 CREATE TABLE YES_DEV.T11 (ID BINARY_BIGINT NOT NULL,NAME VARCH…

【管理咨询宝藏93】大型制造集团数字化转型设计方案

【管理咨询宝藏93】大型制造集团数字化转型设计方案 【格式】PDF版本 【关键词】国际咨询公司、制造型企业转型、数字化转型 【核心观点】 - 235页大型制造型集团数字化转型方案设计&#xff01;细节非常详尽&#xff0c;图表丰富&#xff01; - 系统架构必须采用成熟、具有国…

美食推荐网站设计

**中文摘要&#xff1a;**在当今信息化、网络化的时代背景下&#xff0c;美食文化正逐渐融入人们的日常生活&#xff0c;而网络平台成为人们获取美食信息、分享美食体验的重要途径。为了满足广大美食爱好者对美食信息的探索和推荐需求&#xff0c;本文提出了一种创新的美食推荐…

YOLOv8网络结构介绍

将按照YOLOv8目标检测任务、实例分割任务、关键点检测任务以及旋转目标检测任务的顺序来介绍&#xff0c;主要内容也是在目标检测任务中介绍&#xff0c;其他任务也只是Head层不相同。 1.YOLOv8_det网络结构 首先&#xff0c;YOLOv8网络分成了三部分&#xff0c;分别是主干网络…

Spark云计算平台Databricks使用,创建workspace和Compute计算集群(Spark集群)

Databricks&#xff0c;是属于 Spark 的商业化公司&#xff0c;由美国加州大学伯克利 AMP 实验室的 Spark 大数据处理系统多位创始人联合创立。Databricks 致力于提供基于 Spark 的云服务&#xff0c;可用于数据集成&#xff0c;数据管道等任务。 1 创建workspace 点击创建wor…

word 毕业论文格式调整

添加页眉页脚 页眉 首先在页面上端页眉区域双击&#xff0c;即可出现“页眉和页脚”设置页面&#xff1a; 页眉左右两端对齐 如果想要页眉页脚左右两端对齐&#xff0c;可以选择添加三栏页眉&#xff0c;然后将中间那一栏删除&#xff0c;即可自动实现左右两端对齐&#x…

Spring Boot集成Ldap快速入门Demo

1.Ldap介绍 LDAP&#xff0c;Lightweight Directory Access Protocol&#xff0c;轻量级目录访问协议. LDAP是一种特殊的服务器&#xff0c;可以存储数据数据的存储是目录形式的&#xff0c;或者可以理解为树状结构&#xff08;一层套一层&#xff09;一般存储关于用户、用户…

吴恩达机器学习笔记:第 9 周-17大规模机器学习(Large Scale Machine Learning)17.3-17.4

目录 第 9 周 17、 大规模机器学习(Large Scale Machine Learning)17.3 小批量梯度下降17.4 随机梯度下降收敛 第 9 周 17、 大规模机器学习(Large Scale Machine Learning) 17.3 小批量梯度下降 小批量梯度下降算法是介于批量梯度下降算法和随机梯度下降算法之间的算法&…

基于Springboot的线上教学平台

基于SpringbootVue的线上教学平台设计与实现 开发语言&#xff1a;Java数据库&#xff1a;MySQL技术&#xff1a;SpringbootMybatis工具&#xff1a;IDEA、Maven、Navicat 系统展示 用户登录 首页 学习资料 交流论坛 试卷列表 公告信息 后台登录 后台首页 学员管理 资料类型…

Junit 测试中如何对异常进行断言

本文对在 Junit 测试中如何对异常进行断言的几种方法进行说明。 使用 Junit 5 如果你使用 Junit 5 的话&#xff0c;你可以直接使用 assertThrows 方法来对异常进行断言。 代码如下&#xff1a; Exception exception assertThrows(NumberFormatException.class, () -> {n…

Universal Thresholdizer:将多种密码学原语门限化

参考文献&#xff1a; [LS90] Lapidot D, Shamir A. Publicly verifiable non-interactive zero-knowledge proofs[C]//Advances in Cryptology-CRYPTO’90: Proceedings 10. Springer Berlin Heidelberg, 1991: 353-365.[Shoup00] Shoup V. Practical threshold signatures[C…

七、 数据出境安全评估申报需要多长时间?

《评估申报指南&#xff08;第二版&#xff09;》未区分数据处理者进行数据出境安全评估线上申报和线下申报整体所需时间。一般情况下&#xff0c;数据出境安全评估的申报时长周期如图所示&#xff1a; 根据《评估申报指南&#xff08;第二版&#xff09;》第二条的规定&#…

Spirng-IOC零碎知识点

Spirng IOC 依赖注入 根据名称注入 <?xml version"1.0" encoding"UTF-8"?> <beansxmlns"http://www.springframework.org/schema/beans"xmlns:xsi"http://www.w3.org/2001/XMLSchema-instance"xmlns:util"http://w…

引入RabbitMQ

前置条件 docker 安装 mq docker run \-e RABBITMQ_DEFAULT_USERdudu \-e RABBITMQ_DEFAULT_PASS123456 \-v mq-plugins:/plugins \--name mq \--hostname mq \-p 15672:15672 \-p 5672:5672 \--network hmall \-d \rabbitmq:3.8-management可能会出现&#xff1a;docker: Er…

【深度学习】【Lora训练0】StabelDiffusion,Lora训练,kohya_ss训练

文章目录 环境数据自动标注kohya_ss BLIP2kohya_ss WD14 后续 资源&#xff1a; &#xff08;1&#xff09;训练ui kohya_ss&#xff1a; https://github.com/bmaltais/kohya_ss &#xff08;2&#xff09;kohya_ss 的docker 其他docker https://github.com/ashleykleynhans…

GraphGPT——图结构数据的新语言模型

在人工智能的浪潮中&#xff0c;图神经网络&#xff08;GNNs&#xff09;已经成为理解和分析图结构数据的强大工具。然而&#xff0c;GNNs在面对未标记数据时&#xff0c;其泛化能力往往受限。为了突破这一局限&#xff0c;研究者们提出了GraphGPT&#xff0c;这是一种为大语言…

AcWing 161:电话列表 ← 字典树(Trie 树)之前缀匹配

【题目来源】https://www.acwing.com/problem/content/163/【题目描述】 给出一个电话列表&#xff0c;如果列表中存在其中一个号码是另一个号码的前缀这一情况&#xff0c;那么就称这个电话列表是不兼容的。 假设电话列表如下&#xff1a;Emergency 911 Alice 97625999 Bob …

2022 亚马逊云科技中国峰会,对话开发者论坛

目录 前言 最近整理资料发现还有一些前 2 年的内容没发出来&#xff0c;故补发记录&#xff0c;每年都有新的感悟。 开发者论坛 1. 你认为什么是开发者社区&#xff0c;如何定义一个成功的开发者社区&#xff1f; 我认为可以把开发者社区看成一个 “产品” 来对待&#xff…