使用pixy计算群体遗传学统计量

1 数据过滤

过滤参数:过滤掉次等位基因频率(minor allele frequency,MAF)低于0.05、哈达-温伯格平衡(Hardy– Weinberg equilibrium,HWE)对应的P值低于1e-10或杂合率(heterozygosity rates)偏差过大(± 3 SD)的位点:
去除杂合率(heterozygosity rates)偏差过大(± 3 SD)的个体:
假设,基于Plink计算结果,需要移除sample1(高杂合或低杂合):

#vcftools version:
nohup vcftools --vcf snps_filtered.vcf --remove-indels --maf 0.05 --hwe 1e-10 --max-missing 0.8 --min-meanDP 20 --max-meanDP 500 --remove-indv sample1 --recode --stdout > snps_maf0_05_hwe1e-10_missing0_8.vcf &

vcftools生成的文件中会包含命令行输出,使用sed移除:

nohup sed -i '1,30d' snps_maf0_05_hwe1e-10_missing0_8.vcf &

压缩:

bgzip snps_maf0_05_hwe1e-10_missing0_8.vcf
tabix snps_maf0_05_hwe1e-10_missing0_8.vcf.gz

2 计算 F S T 、 D X Y 、 P i F_{ST}、D_{XY}、Pi FSTDXYPi

安装软件包


nohup pixy --stats pi fst dxy --vcf snps_maf0_05_hwe1e-10_missing0_8.vcf.gz --populations pop.txt --window_size 10000 --bypass_invariant_check 'yes' --n_cores 15 --output_folder results &

3 可视化

可视化之前需要将染色体编号替换为数值:

bash ~/gaoyue/GWAs/script/chr_tran.sh raw_results/pixy_dxy.txt results/pixy_dxy.txt
bash ~/gaoyue/GWAs/script/chr_tran.sh raw_results/pixy_fst.txt results/pixy_fst.txt
bash ~/gaoyue/GWAs/script/chr_tran.sh raw_results/pixy_pi.txt results/pixy_pi.txt
#load packages:
library(ggplot2)
library(tidyverse)

#---------------------------------------------------------------------------------#
#             1.define a function to convert the pixy outputs                     #
#---------------------------------------------------------------------------------#
pixy_to_long <- function(pixy_files){

  pixy_df <- list()

  for(i in 1:length(pixy_files)){

    stat_file_type <- gsub(".*_|.txt", "", pixy_files[i])

    if(stat_file_type == "pi"){

      df <- read_delim(pixy_files[i], delim = "\t")
      df <- df %>%
        gather(-pop, -window_pos_1, -window_pos_2, -chromosome,
               key = "statistic", value = "value") %>%
        rename(pop1 = pop) %>%
        mutate(pop2 = NA)

      pixy_df[[i]] <- df


    } else{

      df <- read_delim(pixy_files[i], delim = "\t")
      df <- df %>%
        gather(-pop1, -pop2, -window_pos_1, -window_pos_2, -chromosome,
               key = "statistic", value = "value")
      pixy_df[[i]] <- df

    }

  }

  bind_rows(pixy_df) %>%
    arrange(pop1, pop2, chromosome, window_pos_1, statistic)

}

#---------------------------------------------------------------------------------#
#                      2.use new function we just defined:                        #
#---------------------------------------------------------------------------------#
## Rcau则替换为对应的文件夹
pixy_folder <- "/nfs_fs/nfs3/gaoyue/gaoyue/Fst/Rdeb_Fst/results/"
pixy_files <- list.files(pixy_folder, full.names = TRUE)
pixy_df <- pixy_to_long(pixy_files)

#---------------------------------------------------------------------------------#
#                                      3.plot:                                    #
#---------------------------------------------------------------------------------#
# create a custom labeller for special characters in pi/dxy/fst
pixy_labeller <- as_labeller(c(avg_pi = "pi",
                             avg_dxy = "D[XY]",
                             avg_wc_fst = "F[ST]"),
                             default = label_parsed)

# plotting summary statistics across all chromosomes
pixy_df %>%
  mutate(chrom_color_group = case_when(as.numeric(chromosome) %% 2 != 0 ~ "even",
                                 chromosome == "X" ~ "even",
                                 TRUE ~ "odd" )) %>%
  mutate(chromosome = factor(chromosome, levels = c(1:22, "X", "Y"))) %>%
  filter(statistic %in% c("avg_pi", "avg_dxy", "avg_wc_fst")) %>%
  ggplot(aes(x = (window_pos_1 + window_pos_2)/2, y = value, color = chrom_color_group))+
  geom_point(size = 0.5, alpha = 0.5, stroke = 0)+
  facet_grid(statistic ~ chromosome,
             scales = "free_y", switch = "x", space = "free_x",
             labeller = labeller(statistic = pixy_labeller,
                                 value = label_value))+
  xlab("Chromsome")+
  ylab("Statistic Value")+
  scale_color_manual(values = c("grey50", "black"))+
  theme_classic()+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.spacing = unit(0.1, "cm"),
        strip.background = element_blank(),
        strip.placement = "outside",
        legend.position ="none")+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,NA))

在这里插入图片描述

Ending!

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

两万字图文详解!InnoDB锁专题!

前言 本文将跟大家聊聊 InnoDB 的锁。本文比较长&#xff0c;包括一条 SQL 是如何加锁的&#xff0c;一些加锁规则、如何分析和解决死锁问题等内容&#xff0c;建议耐心读完&#xff0c;肯定对大家有帮助的。 为什么需要加锁呢&#xff1f; InnoDB 的七种锁介绍 一条 SQL 是…

2023.11.16-hive sql高阶函数lateral view,与行转列,列转行

目录 0.lateral view简介 1.行转列 需求1: 需求2: 2.列转行 解题思路: 0.lateral view简介 hive函数 lateral view 主要功能是将原本汇总在一条&#xff08;行&#xff09;的数据拆分成多条&#xff08;行&#xff09;成虚拟表&#xff0c;再与原表进行笛卡尔积&#xff0c…

那些让我苦笑不得的 Bug:编码之路的坎坷经历

文章目录 1. CSS 中的样式“消失”问题2. JavaScript 的变量命名引发的混乱3. 时间格式的困扰4. 数据库查询条件引发的错误结语 &#x1f389;欢迎来到Java学习路线专栏~那些让我苦笑不得的 Bug&#xff1a;编码之路的坎坷经历 ☆* o(≧▽≦)o *☆嗨~我是IT陈寒&#x1f379;✨…

了解一下知识付费系统的开发流程和关键技术点

知识付费系统的开发既涉及到前端用户体验&#xff0c;又需要强大的后端支持和复杂的付费逻辑。在这篇文章中&#xff0c;我们将深入探讨知识付费系统的开发流程和关键技术点&#xff0c;并提供一些相关的技术代码示例。 1. 需求分析和规划&#xff1a; 在着手开发知识付费系…

小学生加减乘除闯关运算练习流量主微信小程序开发

小学生加减乘除闯关运算练习流量主微信小程序开发 经过本次更新&#xff0c;我们增加了新的功能和特性&#xff0c;以提升用户体验和运算练习的趣味性&#xff1a; 能量石与激励视频&#xff1a;用户可以通过观看激励视频来获取能量石&#xff0c;这些能量石可以用于解锁收费…

ctfshow 文件上传 151-161

文件上传也好久没做了。。 手很生了 151 前端绕过 只能上传png文件 使用bp抓包&#xff0c;修改文件名后缀为php 上传成功&#xff0c;发现文件上传路径 使用蚁剑连接 找到flag 152 152 后端校验 跟上一关一样 表示后面即使执行错误&#xff0c;也不报错 抓包修改文件…

招投标系统软件源码,招投标全流程在线化管理

功能描述 1、门户管理&#xff1a;所有用户可在门户页面查看所有的公告信息及相关的通知信息。主要板块包含&#xff1a;招标公告、非招标公告、系统通知、政策法规。 2、立项管理&#xff1a;企业用户可对需要采购的项目进行立项申请&#xff0c;并提交审批&#xff0c;查看所…

TypeError: Cannot read properties of undefined (reading ‘0‘)

1、在使用<el-dropdown>会报这个错误 原因&#xff1a;使用v-if控制显隐&#xff0c;找不到该节点就会开始报错 解决&#xff1a;使用v-show就可以了

自定义控件封装

上边对两个控件整合时用函数指针是因为QSpinBox::valueChange有重载版本 自定义的接口放在类外 你设计的界面可以通过提升为调用这些接口 添加qt设计师界面类

winform窗体、控件的简单封装,重做标题栏

1标题栏封装中使用了以下技术&#xff1a; 知识点&#xff1a; 1.父类、子类的继承&#xff1b; 2.窗体之间的继承&#xff1b; 3.自定义控件的绘制&#xff1b; 4.多线程在窗体间的应用&#xff1b; 5.窗体和控件的封装&#xff1b; 6.回调函数&#xff08;委托&#xff09;&…

Redis快速入门(基础篇)

简介&#xff1a; 是一个高性能的 key-value数据库。 存在内存中 与其他 key-value 缓存产品有以下三个特点&#xff1a; Redis支持数据的持久化&#xff0c;可以将内存中的数据保持在磁盘中&#xff0c;重启的时候可以再次加载进行使用。 Redis不仅仅支持简单的key-value类…

自动备份pgsql数据库

bat文件中的内容&#xff1a; PATH D:\Program Files\PostgreSQL\13\bin;D:\Program Files\7-Zip set PGPASSWORD**** pg_dump -h 8.134.151.187 -p 5466 -U sky -d mip_db --schema-only -f D:\DB\backup\%TODAY%-schema-mip_db_ali.sql pg_dump -h 8.134.151.187 -p 5466…

BUUCTF 面具下的flag 1

BUUCTF:https://buuoj.cn/challenges 题目描述&#xff1a; 下载附件&#xff0c;得到一张.jpg图片。 密文&#xff1a; 解题思路&#xff1a; 1、将图片放到Kali中&#xff0c;使用binwalk检测出隐藏zip包。 使用foremost提取zip压缩包到output目录下 解压zip压缩包&…

python自动化第一篇—— 带图文的execl的自动化合并

简述 最近接到一个需求&#xff0c;需要为公司里的一个部门提供一个文件上传自动化合并的系统&#xff0c;以供用户稽核&#xff0c;谈到自动化&#xff0c;肯定是选择python&#xff0c;毕竟python的轮子多。比较了市面上几个用得多的python库&#xff0c;我最终选择了xlwings…

【自然语言处理(NLP)实战】LSTM网络实现中文文本情感分析(手把手与教学超详细)

目录 引言&#xff1a; 1.所有文件展示&#xff1a; 1.中文停用词数据&#xff08;hit_stopwords.txt)来源于&#xff1a; 2.其中data数据集为chinese_text_cnn-master.zip提取出的文件。点击链接进入github&#xff0c;点击Code、Download ZIP即可下载。 2.安装依赖库&am…

聊聊模糊测试,以及几种模糊测试工具的介绍!

以下为作者观点&#xff1a; 在当今的数字环境中&#xff0c;漏洞成为攻击者利用系统漏洞的通道&#xff0c;对网络安全构成重大威胁。这些漏洞可能存在于硬件、软件、协议实施或系统安全策略中&#xff0c;允许未经授权的访问并破坏系统的完整性。 根据 "常见漏洞与暴露…

在Linux中nacos集群模式部署

一、安装 配置nacos 在Linux中建立一个nacos文件夹 mkdir nacos 把下载的压缩包拉入刚才创建好的nacos文件中 解压 tar -zxvf nacos-server-1.4.1\.tar.gz 修改配置文件 进入nacos文件中的conf文件的cluster.conf.example 修改cluster.conf.example文件 vim cluster.conf.exa…

Django(六、模板层)

文章目录 模板传值模板语法传值特性 模板语法之过滤器常用的过滤器模板层之标签模板中的标签的格式为标签之if判断 标签之for循环模板的继承与导入模板导入导入格式 模板传值 """ 模板层三种语法 {{}}:主要与数据值相关 {%%}:主要与逻辑相关 {##}&#xff1a;模…