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

2024-05-10 07:52

本文主要是介绍生信分析进阶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)

矫正后结果

矫正后文件

这篇关于生信分析进阶2 - 利用GC含量的Loess回归矫正reads数量的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/975878

相关文章

Spring Security 从入门到进阶系列教程

Spring Security 入门系列 《保护 Web 应用的安全》 《Spring-Security-入门(一):登录与退出》 《Spring-Security-入门(二):基于数据库验证》 《Spring-Security-入门(三):密码加密》 《Spring-Security-入门(四):自定义-Filter》 《Spring-Security-入门(五):在 Sprin

Java进阶13讲__第12讲_1/2

多线程、线程池 1.  线程概念 1.1  什么是线程 1.2  线程的好处 2.   创建线程的三种方式 注意事项 2.1  继承Thread类 2.1.1 认识  2.1.2  编码实现  package cn.hdc.oop10.Thread;import org.slf4j.Logger;import org.slf4j.LoggerFactory

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题,本文将一起从代码层去分析为什么没有建立索引? 开源ERP项目地址:https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置

C#实战|大乐透选号器[6]:实现实时显示已选择的红蓝球数量

哈喽,你好啊,我是雷工。 关于大乐透选号器在前面已经记录了5篇笔记,这是第6篇; 接下来实现实时显示当前选中红球数量,蓝球数量; 以下为练习笔记。 01 效果演示 当选择和取消选择红球或蓝球时,在对应的位置显示实时已选择的红球、蓝球的数量; 02 标签名称 分别设置Label标签名称为:lblRedCount、lblBlueCount

[MySQL表的增删改查-进阶]

🌈个人主页:努力学编程’ ⛅个人推荐: c语言从初阶到进阶 JavaEE详解 数据结构 ⚡学好数据结构,刷题刻不容缓:点击一起刷题 🌙心灵鸡汤:总有人要赢,为什么不能是我呢 💻💻💻数据库约束 🔭🔭🔭约束类型 not null: 指示某列不能存储 NULL 值unique: 保证某列的每行必须有唯一的值default: 规定没有给列赋值时的默认值.primary key:

SWAP作物生长模型安装教程、数据制备、敏感性分析、气候变化影响、R模型敏感性分析与贝叶斯优化、Fortran源代码分析、气候数据降尺度与变化影响分析

查看原文>>>全流程SWAP农业模型数据制备、敏感性分析及气候变化影响实践技术应用 SWAP模型是由荷兰瓦赫宁根大学开发的先进农作物模型,它综合考虑了土壤-水分-大气以及植被间的相互作用;是一种描述作物生长过程的一种机理性作物生长模型。它不但运用Richard方程,使其能够精确的模拟土壤中水分的运动,而且耦合了WOFOST作物模型使作物的生长描述更为科学。 本文让更多的科研人员和农业工作者

MOLE 2.5 分析分子通道和孔隙

软件介绍 生物大分子通道和孔隙在生物学中发挥着重要作用,例如在分子识别和酶底物特异性方面。 我们介绍了一种名为 MOLE 2.5 的高级软件工具,该工具旨在分析分子通道和孔隙。 与其他可用软件工具的基准测试表明,MOLE 2.5 相比更快、更强大、功能更丰富。作为一项新功能,MOLE 2.5 可以估算已识别通道的物理化学性质。 软件下载 https://pan.quark.cn/s/57

【Linux 从基础到进阶】Ansible自动化运维工具使用

Ansible自动化运维工具使用 Ansible 是一款开源的自动化运维工具,采用无代理架构(agentless),基于 SSH 连接进行管理,具有简单易用、灵活强大、可扩展性高等特点。它广泛用于服务器管理、应用部署、配置管理等任务。本文将介绍 Ansible 的安装、基本使用方法及一些实际运维场景中的应用,旨在帮助运维人员快速上手并熟练运用 Ansible。 1. Ansible的核心概念

Flutter 进阶:绘制加载动画

绘制加载动画:由小圆组成的大圆 1. 定义 LoadingScreen 类2. 实现 _LoadingScreenState 类3. 定义 LoadingPainter 类4. 总结 实现加载动画 我们需要定义两个类:LoadingScreen 和 LoadingPainter。LoadingScreen 负责控制动画的状态,而 LoadingPainter 则负责绘制动画。

【编程底层思考】垃圾收集机制,GC算法,垃圾收集器类型概述

Java的垃圾收集(Garbage Collection,GC)机制是Java语言的一大特色,它负责自动管理内存的回收,释放不再使用的对象所占用的内存。以下是对Java垃圾收集机制的详细介绍: 一、垃圾收集机制概述: 对象存活判断:垃圾收集器定期检查堆内存中的对象,判断哪些对象是“垃圾”,即不再被任何引用链直接或间接引用的对象。内存回收:将判断为垃圾的对象占用的内存进行回收,以便重新使用。