重测序数据分析流程丨操作步骤与代码与代码脚本
群体重测序数据分析笔记
在生物信息学中,群体重测序数据的挖掘和分析对于理解生物的进化、自然选择以及功能基因的定位等研究具有重要的意义。今天分享的笔记内容是群体遗传学相关的知识点,下面将一步步介绍整个重测序分析的流程和方法。
分析常用流程和方法简述
常用的重测序分析流程一般包含以下步骤:
-
质控和数据准备:
这一步包括对原始测序数据进行质量评估和控制,以及数据格式的转换,常用流程包括测序数据的质控(QC)、比对、标记重复、排序、建立索引和变异检测。常用的软件工具包括FastQC、BWA、SAMtools、Picard和GATK等。
首先,我们需要进行质量控制,确保我们的数据是准确可靠的。常用的质量控制工具有FastQC和Trimmomatic。以下是用Trimmomatic进行质量控制的例子:
java -jar trimmomatic.jar PE -phred33 \
input_forward.fq.gz input_reverse.fq.gz \
output_forward_paired.fq.gz output_forward_unpaired.fq.gz \
output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:15 MINLEN:36
-
序列比对:
将处理好的数据比对到参考基因组上,得到比对结果,比对可以使用BWA这个工具。以下是比对命令的示例:
# 建立参考基因组索引
bwa index ref.fa
# 比对
bwa mem ref.fa read1.fq read2.fq > aln.sam
-
变异检测:
根据比对结果,检测和分析基因序列中的变异,变异检测通常使用SAMtools和BCFtools。以下是如何使用这些工具进行变异检测的示例:
# 转换格式
samtools view -S -b aln.sam > aln.bam
# 排序
samtools sort aln.bam -o aln_sorted.bam
# 检测变异
samtools mpileup -uf ref.fa aln_sorted.bam | bcftools call -cv - > var.raw.vcf
-
变异注释和分析:
对检测到的变异进行注释和深度分析,包括群体结构分析、选择性消除分析、全基因组关联分析等。
以下为流程简单示意
# 质量控制
fastqc raw_data.fq
# 比对
bwa mem reference.fa raw_data.fq > aln.sam
# 标记重复
picard MarkDuplicates I=aln.sam O=marked.bam M=metrics.txt
# 排序
samtools sort -O bam -o sorted.bam -T temp.prefix marked.bam
# 建立索引
samtools index sorted.bam
# 变异检测
gatk HaplotypeCaller -R reference.fa -I sorted.bam -O variants.vcf
上游分析变异检测方法代码
在进行测序数据序列文件上游变异检测时,我们通常使用如GATK工具。以下是一段使用GATK进行SNP和Indel检测的代码:
# SNP检测
gatk --java-options "-Xmx4g" \
HaplotypeCaller -R reference.fa \
-I sorted.bam \
-O raw_snps.vcf
# Indel检测
gatk --java-options "-Xmx4g" \
HaplotypeCaller -R reference.fa \
-I sorted.bam --ploidy 2 \
--genotyping-mode DISCOVERY \
-stand_emit_conf 10 \
-stand_call_conf 30 \
-O raw_indels.vcf
群体结构分析方法
对于群体结构的分析,我们通常会使用程序如PLINK和ADMIXTURE。PLINK可以用于生成适合ADMIXTURE分析的数据,而ADMIXTURE可以进行群体结构分析。
# PLINK生成适用于ADMIXTURE的数据
plink --file input_data \
--make-bed --out plink_output
# ADMIXTURE进行群体结构分析
admixture --cv plink_output.bed K
变异位点选择性消除分析
在基因组变异位点选择性消除分析中,我们可以使用vcftools和R等工具来评估基因位点多态性和选择分化差异。以下是具体的操作:
# 使用vcftools计算窗口内的多态性
vcftools --vcf var.raw.vcf --window-pi 50000
# 在R中进行选择分化差异判断
# 以下是示例代码,具体代码需要根据实际情况编写
library(ape)
data <- read.table("fst.txt", header = T)
fst <- data$V4
sig_level <- qnorm(1 - 0.05 / 2) / sqrt(2)
outliers <- which(fst > sig_level)
GWAS全基因组关联分析
GWAS全基因组关联分析我们可以使用GAPIT包进行分析。以下是用R语言进行GWAS全基因组关联分析的示例代码:
library(GAPIT)
genotype_file <- "genotype.hmp.txt"
phenotype_file <- "phenotype.txt"
GAPIT_data <- GAPIT.Data(fileHapmap = genotype_file,
filePhenotype = phenotype_file)
GAPIT_GLM <- GAPIT.GLM(GAPIT_data)
GAPIT_MLM <- GAPIT.MLM(GAPIT_data)
PCA分析与进化树绘制
我们可以使用plink和gcta进行PCA分析,并使用ggtree进行群体结构进化树的绘制。以下是具体的操作:
# 使用plink和gcta进行PCA分析
plink --bfile plink --make-grm-bin --out plink
gcta --grm-bin plink --pca 10 --out plink
# 在R中使用ggtree绘制进化树
# 以下是示例代码,具体代码需要根据实际情况编写
library(ggtree)
tree <- read.tree("tree.nwk")
ggtree(tree) + geom_tiplab()
从vcf文件中挖掘显著变异位点的频率变化信息
我们可以使用vcftools工具从vcf文件中挖掘显著变异位点的频率变化信息,以下是具体的操作:
vcftools --vcf var.raw.vcf --freq --out freq
以上就是群体重测序数据的挖掘与分析的整个过程,每一步都需要我们仔细和认真的处理。
总结
以下是一个简单的bash脚本,实现对多个样品的批量重测序分析,这个脚本使用了一个循环来处理每一个样品。请注意这只是一个示例脚本,您可能需要根据实际情况对其进行修改或优化。
#!/bin/bash
# 路径参数
REF_GENOME_PATH=/path/to/your/reference/genome
SAMPLE_LIST=/path/to/your/sample/list
SAMPLE_PATH=/path/to/your/sample/data
WORKING_DIR=/path/to/your/working/directory
# 创建工作目录
mkdir -p ${WORKING_DIR}
# 循环处理每一个样品
while read SAMPLE; do
echo "Processing sample ${SAMPLE}..."
# 1. 质控和数据准备
java -jar trimmomatic.jar PE -phred33 \
${SAMPLE_PATH}/${SAMPLE}_1.fq.gz \
${SAMPLE_PATH}/${SAMPLE}_2.fq.gz \
${WORKING_DIR}/${SAMPLE}_1_paired.fq.gz \
${WORKING_DIR}/${SAMPLE}_1_unpaired.fq.gz \
${WORKING_DIR}/${SAMPLE}_2_paired.fq.gz \
${WORKING_DIR}/${SAMPLE}_2_unpaired.fq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
# 2. 序列比对
bwa mem ${REF_GENOME_PATH} \
${WORKING_DIR}/${SAMPLE}_1_paired.fq.gz \
${WORKING_DIR}/${SAMPLE}_2_paired.fq.gz > ${WORKING_DIR}/${SAMPLE}.sam
# 3. 变异检测
samtools view -S -b ${WORKING_DIR}/${SAMPLE}.sam > ${WORKING_DIR}/${SAMPLE}.bam
samtools sort ${WORKING_DIR}/${SAMPLE}.bam -o ${WORKING_DIR}/${SAMPLE}_sorted.bam
samtools mpileup -uf ${REF_GENOME_PATH} ${WORKING_DIR}/${SAMPLE}_sorted.bam | bcftools call -cv - > ${WORKING_DIR}/${SAMPLE}_var.raw.vcf
# 4. 提取VCF文件中的频率信息
vcftools --vcf ${WORKING_DIR}/${SAMPLE}_var.raw.vcf --freq --out ${WORKING_DIR}/${SAMPLE}_freq
echo "Finished processing sample ${SAMPLE}."
done < ${SAMPLE_LIST}
在这个脚本中,我们假设有一个文本文件 ${SAMPLE_LIST} 包含所有待处理的样品名称,每个样品对应一对fastq.gz文件,文件名分别为 ${SAMPLE}_1.fq.gz 和 ${SAMPLE}_2.fq.gz。所有这些文件都存储在 ${SAMPLE_PATH} 路径下。处理过程中的中间文件和结果文件都会存储在 ${WORKING_DIR} 路径下。
使用这个脚本之前,你需要修改上面的四个路径参数,使它们指向实际的路径。另外,这个脚本只执行了一部分的分析步骤,如果你需要执行更多的步骤,你可以在脚本中添加相应的命令。
本文由 mdnice 多平台发布
相关文章:
重测序数据分析流程丨操作步骤与代码与代码脚本
群体重测序数据分析笔记 在生物信息学中,群体重测序数据的挖掘和分析对于理解生物的进化、自然选择以及功能基因的定位等研究具有重要的意义。今天分享的笔记内容是群体遗传学相关的知识点,下面将一步步介绍整个重测序分析的流程和方法。 分析常用流程和…...
npm -v无法显示版本号
情况: 删除C盘下.npmrc文件后解决。路径 C:\Users\Dell 记录一下这个解法。...
【Vue】父子组件值及方法传递使用
父子组件值、方法引用 1、值 1.1 父组件获取子组件值 父组件 <template><div><button click"getChildValue">click</button><child ref"child"></child></div> </template><script> import Child…...
医药化工企业洁净厂房改造消防防爆安全的重要性
设计 【摘要】:近年来,我国医药化工企业规模不断扩大。医药化工企业的情况复杂,稍有不慎将发生火灾或者爆炸,对人员生命以及财产安全造成巨大的损害,酿成悲剧。所以,“三同时”原则的落实,如何…...
Web开发中防止SQL注入
一、SQL注入简介 SQL注入是比较常见的网络攻击方式之一,它不是利用操作系统的BUG来实现攻击,而是针对程序员编写时的疏忽,通过SQL语句,实现无账号登录,甚至篡改数据库。 二、SQL注入攻击的总体思路 1.寻找到SQL注入…...
【LeetCode-中等】剑指 Offer 35. 复杂链表的复制(详解)
目录 题目 方法1:错误的方法(初尝试) 方法2:复制、拆开 方法3:哈希表 总结 题目 请实现 copyRandomList 函数,复制一个复杂链表。在复杂链表中,每个节点除了有一个 next 指针指向下一个节…...
QT图形视图系统 - 使用一个项目来学习QT的图形视图框架 -第一篇
文章目录 QT图形视图系统介绍开始搭建MainWindow框架设置scene的属性缩放功能的添加加上标尺 QT图形视图系统 介绍 详细的介绍可以看QT的官方助手,那里面介绍的详细且明白,需要一定的英语基础,我这里直接使用一个开源项目来介绍QGraphicsVi…...
Cat.1如何成为物联网业务加速器?
随着Cat.1芯片及模组在功耗和成本上的不断优化,在窄带物联网领域,越来越多的终端客户把Cat.1当做与NB-IoT相比较的第二选择。越来越多的表计、烟感、市政等行业终端将Cat.1模组应用于非集中化部署的上报类终端业务中,Cat.1这只“网红猫”仍保…...
Qt应用开发(基础篇)——布局管理 Layout Management
目录 一、前言 二:相关类 三、水平、垂直、网格和表单布局 四、尺寸策略 一、前言 在实际项目开发中,经常需要使用到布局,让控件自动排列,不仅节省控件还易于管控。Qt布局系统提供了一种简单而强大的方式来自动布局小部件中的…...
Python web实战之 Django 的 ORM 框架详解
本文关键词:Python、Django、ORM。 概要 在 Python Web 开发中,ORM(Object-Relational Mapping,对象关系映射)是一个非常重要的概念。ORM 框架可以让我们不用编写 SQL 语句,就能够使用对象的方式来操作数据…...
pycharm制作柱状图
Bar - Bar_rotate_xaxis_label 解决标签名字过长的问题 from pyecharts import options as opts from pyecharts.charts import Barc (Bar().add_xaxis(["高等数学1,2","C语言程序设计","python程序设计","大数据导论",…...
静态资源导入探究
静态资源可以在哪里找呢?我们看看源码 从这个类进去 里面有个静态类 WebMvcAutoConfigurationAdapter 有个配置类,将这个类的对象创建并导入IOC容器里 这个静态类下有个方法 addResourceHandlers(ResourceHandlerRegistry registry)静态资源处理器 若自…...
安全狗V3.512048版本绕过
安全狗安装 安全狗详细安装、遇见无此服务器解决、在windows中命令提示符中进入查看指定文件夹手动启动Apache_安全狗只支持 glibc_2.14 但是服务器是2.17_黑色地带(崛起)的博客-CSDN博客 安全狗 safedogwzApacheV3.5.exe 右键电脑右下角安全狗图标-->选择插件-->安装…...
prometheus监控k8s kube-proxy target down
prometheus kube-proxy target down 解决 修改配置 kubectl edit cm/kube-proxy -n kube-systemmetricsBindAddress: "0.0.0.0:10249"删除 kube-proxy pod 使之重启应用配置 kubectl delete pod --force `kubectl get pod -n kube-system |grep kube-proxy|awk {pr…...
SPSS数据分析--假设检验的两种原假设取舍决定方式
假设检验的两种原假设取舍决定方式 在t检验,相关分析,回归分析,方差分析,卡方检验等等分析方法中,都需要用到假设检验。假设检验的步骤一般如下: 提出假设:H0 vs H1 ;假设原假设H0 成立的情况…...
Python实现猫狗分类
不废话了,直接上代码: def load_imagepath_from_csv(csv_name):image_path []with open(csv_name,r) as file:csv_reader csv.reader(file)next(csv_reader)for row in csv_reader:image_path.append(row[0])return image_pathimport csv csv_name &…...
pjsip、pjsua2+bcg729 windows下编译java版本
文章目录 简要说明流程步骤 简要说明 基本参考的这里 https://docs.pjsip.org/en/latest/get-started/windows/build_instructions.html#building-the-projects 我这里主要是为了生成pjsua2.dll 用于在java下调用。 其中 libbcg729.dll 是通过vcpkg来进行安装。 pjsip使用vs2…...
尝试多数据表 sqlite
C 唯一值得骄傲的地方就是 通过指针来回寻址 😂 提高使用的灵活性 小脚本buff 加成...
Keil出现Flash Timeout.Reset the Target and try it again.我有一种解决方法
2.解决方法 网上查找了找原因,是因为之前代码设置了读保护功能。 读保护即大家通常说的“加密”,是作用于整个Flash存储区域。一旦设置了Flash的读保护,内置的Flash存储区只能通过程序的正常执行才能读出,而不能通过下述任何一种…...
纯粹即刻,畅享音乐搜索的轻松体验
纯粹即刻,畅享音乐搜索的轻松体验 在当今快节奏的生活中,我们常常渴望一种简单而便捷的方式来探索和享受音乐。现在,你可以纯粹即刻地畅享音乐搜索的轻松体验。无论你是寻找热门歌曲还是探索不同风格的音乐,这款应用将为你带来随…...
KubeSphere 容器平台高可用:环境搭建与可视化操作指南
Linux_k8s篇 欢迎来到Linux的世界,看笔记好好学多敲多打,每个人都是大神! 题目:KubeSphere 容器平台高可用:环境搭建与可视化操作指南 版本号: 1.0,0 作者: 老王要学习 日期: 2025.06.05 适用环境: Ubuntu22 文档说…...
网络编程(Modbus进阶)
思维导图 Modbus RTU(先学一点理论) 概念 Modbus RTU 是工业自动化领域 最广泛应用的串行通信协议,由 Modicon 公司(现施耐德电气)于 1979 年推出。它以 高效率、强健性、易实现的特点成为工业控制系统的通信标准。 包…...
label-studio的使用教程(导入本地路径)
文章目录 1. 准备环境2. 脚本启动2.1 Windows2.2 Linux 3. 安装label-studio机器学习后端3.1 pip安装(推荐)3.2 GitHub仓库安装 4. 后端配置4.1 yolo环境4.2 引入后端模型4.3 修改脚本4.4 启动后端 5. 标注工程5.1 创建工程5.2 配置图片路径5.3 配置工程类型标签5.4 配置模型5.…...
相机Camera日志实例分析之二:相机Camx【专业模式开启直方图拍照】单帧流程日志详解
【关注我,后续持续新增专题博文,谢谢!!!】 上一篇我们讲了: 这一篇我们开始讲: 目录 一、场景操作步骤 二、日志基础关键字分级如下 三、场景日志如下: 一、场景操作步骤 操作步…...
centos 7 部署awstats 网站访问检测
一、基础环境准备(两种安装方式都要做) bash # 安装必要依赖 yum install -y httpd perl mod_perl perl-Time-HiRes perl-DateTime systemctl enable httpd # 设置 Apache 开机自启 systemctl start httpd # 启动 Apache二、安装 AWStats࿰…...
java调用dll出现unsatisfiedLinkError以及JNA和JNI的区别
UnsatisfiedLinkError 在对接硬件设备中,我们会遇到使用 java 调用 dll文件 的情况,此时大概率出现UnsatisfiedLinkError链接错误,原因可能有如下几种 类名错误包名错误方法名参数错误使用 JNI 协议调用,结果 dll 未实现 JNI 协…...
2021-03-15 iview一些问题
1.iview 在使用tree组件时,发现没有set类的方法,只有get,那么要改变tree值,只能遍历treeData,递归修改treeData的checked,发现无法更改,原因在于check模式下,子元素的勾选状态跟父节…...
Linux 内存管理实战精讲:核心原理与面试常考点全解析
Linux 内存管理实战精讲:核心原理与面试常考点全解析 Linux 内核内存管理是系统设计中最复杂但也最核心的模块之一。它不仅支撑着虚拟内存机制、物理内存分配、进程隔离与资源复用,还直接决定系统运行的性能与稳定性。无论你是嵌入式开发者、内核调试工…...
毫米波雷达基础理论(3D+4D)
3D、4D毫米波雷达基础知识及厂商选型 PreView : https://mp.weixin.qq.com/s/bQkju4r6med7I3TBGJI_bQ 1. FMCW毫米波雷达基础知识 主要参考博文: 一文入门汽车毫米波雷达基本原理 :https://mp.weixin.qq.com/s/_EN7A5lKcz2Eh8dLnjE19w 毫米波雷达基础…...
Golang——7、包与接口详解
包与接口详解 1、Golang包详解1.1、Golang中包的定义和介绍1.2、Golang包管理工具go mod1.3、Golang中自定义包1.4、Golang中使用第三包1.5、init函数 2、接口详解2.1、接口的定义2.2、空接口2.3、类型断言2.4、结构体值接收者和指针接收者实现接口的区别2.5、一个结构体实现多…...
