计算fst、pi等群体遗传特征数据
先得到NLR所有基因的vcf信息
将每个基因的信息提取成单独的vcf文件
vcftools进行计算
目 录
一、处理vcf文件
文件位置：
（一）提取信息的代码
1. python循环嵌套注意事项：
2. python中把文件按行读取的方法：
（二）提交命令的代码
1. shell中的列表循环
2. shell中批量提交py脚本
3. Python extend()方法添加元素
4. python writelines(string)写入文件：
5.Linux nohup命令详解
6.shell 截取文件第一列并去重命令
一、处理vcf文件
分别获得Morex的NLR家族基因在snp.vcf和indel.vcf的信息
！文件信息：
文件位置：
/data/database/barley_reseq
输入文件：
morex_info2.csv Hv.final.sum.snp.foreig.eff.chr1H.vcf(共14个)
输出文件：
morex_snp.txt morex_indel.txt gene_id.vcf
（一）提取信息的代码
#判断snp的位置是否在morex的位置信息中
import sys
import pandas as pd
#snp=open(sys.argv[1])
#snp=open("Hv.final.sum.snp.foreig.eff.chr1H.vcf","r")
location={}
Chr={}
Start={}
End={}
gene_start={}
gene_end={}
Info={}
info_list=[]morex=open("morex_info2.csv","r")
for line in morex.readlines():line=line.replace("\n","")chr1=line.split(",")[1]id=line.split(",")[0]start=line.split(",")[2]end=line.split(",")[3]# for k in location.keys():# if int(k) >= int(start) and int(k) <= int(end):Start[id]=int(start)End[id]=int(end)Chr[id]=chr1
morex.close()for k in Chr.keys():info_list=[]with open(sys.argv[1]) as snp:for l in snp:#l=l.split("\t",2)chr2=l.split("\t",2)[0]loc=int(l.split("\t",2)[1])#info=l[2]if Chr[k] == chr2 and Start[k] <= loc <= End[k]:location[k]=locinfo_list.append(l)#gene_end[k]=End[k]Info[k] = info_list#print(len(Info['HORVU.MOREX.r2.1HG0000230']))
output=open("morex_snp.txt","a+")
for i in location.keys():for j in Info[i]:output.write(i+"\t"+j)
output.close()
1. python循环嵌套注意事项：
在进行嵌套循环时，内循环的发生条件必须在循环内部！！！
python for循环嵌套时，内层循环不执行的问题 - 简书 (jianshu.com)https://www.jianshu.com/p/6d6f7ffa1977
当输出字典为文件时，字典的值为列表 ，可以使用双循环，先循环key值，再对键所对应的值进行循环，然后输出
2. python中把文件按行读取的方法：
5种Python逐行读取文件的方式_Python热爱者的博客-CSDN博客_python 逐行读取https://blog.csdn.net/qdPython/article/details/106160272
（二）提交命令的代码
#!/bin/bash
adress="/data/database/barley_reseq"
#list=(Hv2.final2.sum.indel.eff.chr1H.vcf Hv2.final2.sum.indel.eff.chr2H.vcf Hv2.final2.sum.indel.eff.chr3H.vcf Hv2.final2.sum.indel.eff.chr4H.vcf Hv2.final2.sum.indel.eff.chr5H.vcf Hv2.final2.sum.indel.eff.chr6H.vcf Hv2.final2.sum.indel.eff.chr7H.vcf)
list=(Hv.final.sum.snp.foreig.eff.chr1H.vcf Hv.final.sum.snp.foreig.eff.chr2H.vcf Hv.final.sum.snp.foreig.eff.chr3H.vcf Hv.final.sum.snp.foreig.eff.chr4H.vcf Hv.final.sum.snp.foreig.eff.chr5H.vcf Hv.final.sum.snp.foreig.eff.chr6H.vcf Hv.final.sum.snp.foreig.eff.chr7H.vcf)
for i in ${list[*]}
do
python snp.py $adress/$i
wait
done1. shell中的列表循环
在shell脚本中，循环列表中所有元素用list[@] 或 list[*]
shell之列表的定义与循环 - 定静沉行 - 博客园 (cnblogs.com)https://www.cnblogs.com/zyy98877/p/10234527.html
2. shell中批量提交py脚本
批量提交脚本时，可以在循环内部写入wait，表示等上一个变量执行结束再进行下一次循环
Shell批量、顺序执行py脚本_追枫萨的博客-CSDN博客_shell批量运行pythonhttps://blog.csdn.net/m0_38052384/article/details/105199879
3. Python extend()方法添加元素
当然，如果希望不将被追加的列表或元组当成一个整体，而是只追加列表中的元素，则可使用列表提供的 extend() 方法。
b_list = ['a', 30]# 追加元组中的所有元素b_list.extend((-2, 3.1))print(b_list)########结果如下
['a', 30, -2, 3.1]b_list.extend(['C', 'R', 'A'])print(b_list)#########结果如下
['a', 30, -2, 3.1, 'C', 'R', 'A']python怎么向列表中添加多个元素_品易HTTP的博客-CSDN博客_python列表怎么添加多个元素https://blog.csdn.net/m0_51713294/article/details/112389296
4. python writelines(string)写入文件：
weitelines可以将一个列表中的内容读入到文件中
>>>>fobj = open('x','w')
>>>>msg = ['write date\n','to x\n','finish\n']
>>>>fobj.writelines(msg)
>>>>fobj.close()
x内容:
write date
to 3.txt
finishpython 写文件write(string), writelines(list) ，读文件 - xushukui - 博客园 (cnblogs.com)https://www.cnblogs.com/nyist-xsk/p/7365786.html
5.Linux nohup命令详解
kill -9 进程号PID
nohup /root/runoob.sh > runoob.log 2>&1 &
Linux nohup 命令 | 菜鸟教程 (runoob.com)https://www.runoob.com/linux/linux-comm-nohup.html
6.shell 截取文件第一列并去重命令
以空格为分隔符截取文件每行的第一个字符串，并用sort排序，再去掉相同的字符串，将结果输出到另一个文件
cat 2.txt | awk -F " " '{print $1}' | sort | uniq > 11.txt
2022.3.30 重新提取morex中的vcf总文件，原因：需要的格式为id+所有内容，分开时只需要把相同基因的所有内容放入一个文件，加上header即可
cat morex_snp |awk -F "\t" '{print $1}' |sort -u|grep "r2.2H"|wc -lcat morex_snp |awk -F "\t" '{print $1}' |sort |uniq|grep "r2.1H"|wc -l在去重时必需先进行排序再去重，或者用sort -u
shell脚本去重的几种方法_小白的进阶的博客-CSDN博客_shell 去重https://blog.csdn.net/laobai1015/article/details/91455406
二、将Morex的总vcf文件分割成单独vcf文件
(一) 分割代码
##将提取到的snp信息按基因提取到不同文件中id_info=[]
id_snp={}i="HORVU.MOREX.r2.7HG0526890"
snp=open("morex_snp.txt","r")
for l in snp.readlines():id = l.split("\t",1)[0]info = l.split("\t",1)[1]id_info.append(info)id_snp[id]=id_infoif i != id:id_info=[]id_info.append(info)id_snp[id]=id_infoi=id
snp.close()
head=open("header","r")
contect=head.readlines()
print(len(id_snp.keys()))
#print(len(set(ID)))
for k in id_snp.keys():for j in id_snp[k]:output=open(k+".vcf","a+")output.writelines(contect)output.write(k+"\t"+j)output.close()1. python 输出文件名含有变量
直接output.write(变量+“.csv”,"w")即可
2. python 对列表去重
直接用set()，结果是随机的
print(len(set(ID)))Python对列表去重的4种方法 - Python基础教程|Python教程|Python入门 - PythonTab中文网https://www.pythontab.com/html/2017/pythonjichu_1204/1194.html
2022.4.5 将单独的vcf文件进行fst均值计算，并整合成result文件，文件内容为：id名+weight_fst+man_fst。同时，对id名进行替换：r2替换为r3。
三、计算单独vcf文件的fst均值
（一）算均值代码及id转换
#计算所有基因的平均和加权fst
#文件后缀“.weir.fst”
# CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
# chr1H 1 1000 1 -1 -1
# chr1H 722001 723000 1 0.0302719 0.0302719
# chr1H 723001 724000 6 0.216685 0.113797
# chr1H 749001 750000 1 0.060739 0.060739#提出来后还需要转换id，所需文件：indel_V2_blastV3.out snp_V2_blastV3.out
#文件地址：D:\张颖\NLR\all_TPM\WGCNA\NLR_WGCNA\hubGene
import os
import glob
import numpy as npos.chdir("/data/user/ZY/ZY/NLR/NLR_base_info/snp/Morex_snp/indel_vcf/result")
# path="/data/user/ZY/ZY/NLR/NLR_base_info/snp/Morex_snp/indel_vcf/result"
# for infile in glob.glob(os.path.join(path, '*.log')):
# os.remove(infile)
files=os.listdir()
####################计算文件fst均值并修改id名##############
def cal_fst(name):f=open(name,"r")invert_id=open("/data/user/ZY/ZY/NLR/NLR_base_info/snp/Morex_snp\indel_V2_blastV3.out","r")W_fst={}M_fst={}weight=[]mean=[]Invert={}Score={}id=name.rsplit(".",3)[0]#m="HORVU.MOREX.r2.1HG0000230"for l in invert_id:r2_id=l.split("\t")[0].rsplit(".",1)[0]r3_id=l.split("\t")[1].rsplit(".",1)[0]if r2_id not in Invert.keys():Invert[r2_id] = r3_idif id in Invert.keys():id=Invert[id]for line in f:if line.startswith("chr"):line=line.split("\t")weight_fst=line[4]mean_fst=line[5]weight.append(float(weight_fst))mean.append(float(mean_fst))W_fst[id]=np.mean(weight)M_fst[id]=np.mean(mean)f.close()return W_fst,M_fst,id# W_fst,M_fst,id=cal_fst("HORVU.MOREX.r2.1HG0000230.windowed.weir.fst")
# print(id)
##########################输出文件###############################
output=open("all_gene_fst_indel.csv","a+")
output.write("ID"+","+"WEIGHTED_FST"+","+"MEAN_FST"+"\n")
for n in files:W_fst,M_fst,id=cal_fst(n)output.write(id+","+str(W_fst[id])+","+str(M_fst[id])+"\n")
output.close()1、列表求均值
import numpy as np
mean1 = np.mean(list1) #平均值
max1 = np.max(list1) #最大值
min1 = np.min(list1) #最小值要保证列表都为数字
2、删除目录下指定后缀文件
import osfrom os import listdir
my_path = 'C:\Python Pool\Test\'
for file_name in listdir(my_path):if file_name.endswith('.txt'):os.remove(my_path + file_name)讲解Python 中删除文件的几种方法-Python教程-PHP中文网https://www.php.cn/python-tutorials-472380.html
3、读取文件夹下的所有文件
import os
path = "D:/Python34/news" #文件夹目录
files= os.listdir(path) #得到文件夹下的所有文件名称
s = []
for file in files: #遍历文件夹if not os.path.isdir(file): #判断是否是文件夹，不是文件夹才打开f = open(path+"/"+file); #打开文件iter_f = iter(f); #创建迭代器str = ""for line in iter_f: #遍历文件，一行行遍历，读取文本str = str + lines.append(str) #每个文件的文本存到list中
print(s) #打印结果4、报错—— cannot perform reduce with flexible type
列表中的元素，并不是数值型，而是string字符串类型，这导致了np.std()函数的报错
TypeError: cannot perform reduce with flexible type 解决方案_leo安静的博客-CSDN博客https://blog.csdn.net/qq_18999357/article/details/80865982
计算fst、pi等群体遗传特征数据相关推荐
- 【群体遗传】Fst（群体间分化指数）
(1)FSTF_{ST}FST是什么?含义是什么? FSTF_{ST}FST,全称为fixation index,是一种用于衡量群体间分化程度的统计检验量(由Wright's F-statisti ...
- 统计遗传学：第三章，群体遗传
3. 群体遗传 大家好,我是飞哥. 前几天推荐了这本书,可以领取pdf和配套数据代码.这里,我将各个章节介绍一下,总结也是学习的过程. 引文部分是原书的谷歌翻译,正文部分是我的理解. 第一部分基础,分 ...
- 重测序群体遗传进化分析之进化树构建
tree 重测序大家都不陌生,它是检测样本基因组变异(SNP,indel,SV,CNV)的主要手段之一,有了这些变异信息,后续可以做很多分析工作,例如: 遗传群体可以进行遗传图谱构建.BSA分析等:大 ...
- 群体遗传，进化分析利器Popgene分享给大家
Popgene分析方法: SSR为显性分子标记,故将实验分析数据作为单倍型数据进行统计,在某一位点上依据条带的有或无分别记为1或0,输入Excel表中形成二元数据矩阵.在此数据基础上进行如下统计分析: ...
- 从幂律分布到特征数据概率分布——12个常用概率分布
在机器学习领域,概率分布对于数据的认识有着非常重要的作用.不管是有效数据还是噪声数据,如果知道了数据的分布,那么在数据建模过程中会得到很大的启示. 首先,如下图所示8个特征数据概率分布情况(已经做归一 ...
- 当 Kubernetes 遇到机密计算，阿里巴巴如何保护容器内数据的安全？
作者 | 贾之光(甲卓) 阿里巴巴高级开发工程师,专注于 Kubernetes 安全沙箱和机密计算领域,主要参与 Incalvare Containers 社区开发. 8 月 26 日,我们发起了第 ...
- 基于流计算 Oceanus（Flink） CDC 做好数据集成场景
作者:黄龙,腾讯 CSIG 高级工程师 数据时代,企业对技术创新和服务水准的要求不断提高,数据已成为企业极其重要的资产.无论是在在企业数据中台的建设,亦或者是打造一站式数据开发和数据治理的PASS平台 ...
- 群体遗传 | haplotype block | HaploBlocker使用
HaploBlocker是分析单条染色体haplotype blocks,这也很容易理解,haplotype block是同一条染色体的某些区域.因此,分析时,需要按染色体或者Scaffold切分VC ...
- pandas使用isna函数和any函数计算返回dataframe中包含缺失值的数据行（rows with missing values in dataframe）
pandas使用isna函数和any函数计算返回dataframe中包含缺失值的数据行(rows with missing values in dataframe) 目录
最新文章
- Smartmail外贸CRMBuild1.0版使用说明书
- java 配置文件加载_Java加载配置文件类
- 如何清理电脑c盘_【电脑】第一期干货：如何正确清理C盘？
- 4固定在底部_有线鼠标之灵魂伴侣，火线竞技4号RGB鼠标线夹
- 最大子序列求和_最大连续子序列和
- Sqlyog的安装与使用
- abab的四字成语_abab式的四字词语
- Bugku-网站被黑
- intel安装mac os
- java字符串在字符中的位置_Java如何获取字符在字符串中的位置
- IDEA修改项目war包名称
- 裁员潮，带给我的思考
- 有的放矢-电气工程师的工作重心
- RAR Extractor - The Unarchiver Pro for mac(解压缩软件)
- 大学开设大数据专业，都安排了哪些课程？
- 计算机勾兑双绝是谁发明,舌尖上的五粮液——记勾兑大师范玉平(图)
- 什么是爱情——碧海青天BBS
- Erebus 0.5 发布，2D 实时角色扮演游戏
- Linux——Makefile文件
- 大数据应用分析解决方案----图书出版