Vst介绍

Vst是通过计算拷贝数的方差来衡量不同群体之间CNV的分化的一个指标，类似于Fst的概念，可以用来鉴定一些高分化的区域。

计算方法如下：

Vpop1是指群体1的copy number的方差；
Vpop2是指群体2的copy number的方差；
Vtotal是全部个体的copy number的方差；
Npop指的对应群体的个体数。

部分文献中，会按照10K的窗口，2K的步长，去计算全基因组上每个窗口的Vst。而一般我们可以用相应的CNV calling的软件检测出每个CNV的位置和拷贝数。因此，我们把这个过程分成两部分：

python脚本计算每个窗口上的绝对拷贝数。
R脚本计算每个窗口的Vst。

计算每个窗口的绝对拷贝数

因为每个检测CNV软件生成结果的格式不大相同，这里把结果中CNV的位置和绝对拷贝数整理成一个TSV格式的文件。

1.文件准备

2.编写脚本

'''
@Author: Guo
@Date: 2020-05-10 10:07:54
@E-mail: willgyw@126.com
@Description:
@LastEditors:
@LastEditTime: 2020-05-13 14:59:18
'''
import sys,os
import logging
import clicklogging.basicConfig(filename=os.path.basename(__file__).replace('.py','.log'),format='%(asctime)s: %(name)s: %(levelname)s: %(message)s',level=logging.DEBUG,filemode='w')
logging.info(f"The command line is:\n\tpython3 {' '.join(sys.argv)}")@click.command()
@click.option('-l', '--lengthfile', help='Input length of each chromosome file', required=True)
@click.option('-t', '--tsv', help='Input copy number tsvfile', required=True)
@click.option('-w', '--windowsize', type=click.INT, help='windowsize', required=True)
@click.option('-s', '--stepsize', type=click.INT, help='stepsize', required=True)
@click.option('-n', '--num', type=click.INT, help='number of samples', required=True)
@click.option('-o', '--output', help='outputfile',default='copy_number_by_window.tsv')def main(lengthfile, tsv, windowsize, stepsize, num, output):step = int(stepsize)win = int(windowsize)samples_num = int(num)fo = open(output, 'w')length_dict = {}with open(lengthfile, 'r') as lengthf:for line in lengthf:lines = line.strip().split()length_dict[lines[0]] = int(lines[1])with open(tsv, 'r') as vcf:for chr in length_dict:for start in range(0, length_dict.get(chr, 0), step):cn_lst = []for num in range(samples_num):cn_lst.append(0)vcf.seek(0)for line in vcf:lines = line.strip().split()c = lines[0]s = lines[1]e = lines[2]if c == chr:if int(start) < int(e) < (int(start) + int(win)):for num in range(samples_num):cn_lst[num] += round(float(lines[(num+3)].split(':')[1]), 2)else:continueelse:continueif sum(cn_lst) > 0:cn_lst_new = [str(round(x, 2)) for x in cn_lst]fo.write('\t'.join([str(chr), str(start), str(int(start) + int(win))]) + '\t' + '\t'.join(cn_lst_new) + '\n')fo.close()if __name__ == '__main__':main()

脚本利用click模块引入命令行参数，需要提供每个染色体的长度、整理好的tsv格式文件、窗口、步长、样本数等参数。

最后会生成类似下图的文件：

窗口位置和每个窗口的拷贝数。

计算每个窗口的Vst

有了上述的结果文件之后，我们就可以利用R脚本计算每个窗口的Vst了。

R的脚本我就不放了，因为不是我写的，尊重原创。python脚本还有很多可以完善的地方，欢迎大家批评指正。

【Python】利用滑动窗口计算全基因组每个窗口上CNV的拷贝数和Vst相关推荐

Python爬虫人工智能大数据全栈视频史上最全合辑教程分享！
Python爬虫人工智能大数据全栈视频史上最全合辑教程分享! 毫无疑问Python是这两年最火的编程语言,不仅容易上手,且在多个行业都可应用.尤其今年人工智能及大数据的发展,Python将会展现更多的 ...
python利用datetime模块计算时间差
python中通过datetime模块可以很方便的计算两个时间的差,datetime的时间差单位可以是天.小时.秒,甚至是微秒,下面我们就来详细看下datetime的强大功能吧今天写了点东西,要计算 ...
【无标题】python利用公式法计算圆周率
# 计算圆周率(公式法) i = 1 j = 1 s = 0 print('******圆周率公式法******') for i in range(1, 100):print(f'循环{i}次,第一次 ...
Python利用meteinfo来计算后向轨迹
## 上面的基本就可以实现了,接下来来尝试算2019年所有天的轨迹 ## 先做2019年的轨迹 # 先获取所有的气象数据于一个列表中 import ostrain_dir=r'H:\meteorolo ...
python 利用datetime 模块计算时间差、当前时间多加一天、一小时、一分钟等操作
1. 常用操作 from datetime import datetimeaa = datetime.now()print(aa) # datetime.datetime(2021, 9, 24, 1 ...
Python关于1、计算一个三位数每位上的数字的三次方之和 2、计算一个四位数每位上的数字的四次方之和
易基因：全基因组CpG密度和DNA甲基化分析方法比较（MeDIP、RRBS和WGBS）| 研究综述
大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因. CpG密度(CpG density)与各种组织中的DNA甲基化相关.基因组按CpG密度分为:CpG岛(CpG island,CGI).C ...
java窗口如何定时关闭_[Java教程]【温故而知新】Javascript窗口效果（全屏显示窗口、定时关闭窗口）_星空网...
[温故而知新]Javascript窗口效果 (全屏显示窗口.定时关闭窗口) 2014-10-30 0 1.全屏显示窗口全屏显示窗口 2.定时关闭窗口定时关闭窗口来源:<HTML.CSS.J ...
flink 处理迟到数据(Trigger、设置水位线延迟时间、允许窗口处理迟到数据、将迟到数据放入侧输出流、代码示例、迟到数据触发窗口计算重复结果处理)
文章目录前言 1.Trigger 2.处理迟到数据 2.1 设置水位线延迟时间 2.2 允许窗口处理迟到数据 2.3 将迟到数据放入侧输出流 3.实操 3.1 代码示例 3.2 中间遇到的异常 3. ...

【Python】利用滑动窗口计算全基因组每个窗口上CNV的拷贝数和Vst

目录