GWAS位點(diǎn)簡單進(jìn)行峰拆分輸出QTNs并基于GFF關(guān)聯(lián)基因

來自rMVP的關(guān)聯(lián)后的數(shù)據(jù);峰拆分輸出QTNs
getposQTNs.py

import os
import sys
from numpy import median
goin =  sys.argv[1]
ranges = int(sys.argv[2])
file = open(goin,"r")
lines = list(file.readlines())
dic = {}
dicout = {}
for i in lines:
    i = i.strip()
    if i.find("SNP") == -1:
        i = i.replace('"','')
        i = i.split(",")
        chrn = i[1] 
        #print(i)
        dicout[i[1]+"-"+str(i[2])] = [i[5],i[7]]
        if chrn not in dic.keys():
            dic[chrn] = []
            dic[chrn].append(int(i[2]))
        else:
            dic[chrn].append(int(i[2]))



#print(dic)
dicoutmid = {}            
for i in dic.keys():
    
    dicoutmid[i] =[]
    listpos = sorted(dic[i])#,reverse=True)
    #print(listpos)
    for j in range(len(listpos)):
        if j == 0:
            gl = []
            gl.append(listpos[j])

        if j + 1 < len(listpos):
            if listpos[j+1] - listpos[j] < ranges:
                gl.append(listpos[j+1])
            else:
                dicoutmid[i].append(gl)
                gl = []
                gl.append(listpos[j+1])
        else:
            dicoutmid[i].append(gl)
#print(dicoutmid)                
for i in dicoutmid.keys():
    for j in dicoutmid[i]:
        if len(j) %2 == 0:
            j = j[:-1]
        pos = median(j)
        effect = dicout[i+"-"+str(int(round(pos)))][0]
        pvalue = dicout[i+"-"+str(int(round(pos)))][1]
        print(i,int(round(pos)),effect,pvalue)

基于GFF關(guān)聯(lián)基因,以水稻RAP注釋為例
(locus.gff文件;https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_2022-09-01.tar.gz

需要處理一下chromosome編號(hào)保證相同

sed -i 's/chr0//g' locus.gff
sed -i 's/chr//g' locus.gff

mergegeneandpoiresult.py

import sys

goin = sys.argv[2] #上面的QTNs輸出
goingff = sys.argv[1] #locus.gff
ranges = int(sys.argv[3]) #關(guān)聯(lián)區(qū)間

filegff = open(goingff,"r")
filel = filegff.readlines()
filegff.close()

dicg ={}
for i in filel:
    i = i.strip()
    i = i.split()
    chrg = i[0]
    gst =  int(i[3])
    ged =  int(i[4])
    chr_locus= chrg +"_"+str(gst) +"_"+ str(ged)
    nameID = i[8].split("ID=")[1].split(";")[0]
    dicg[nameID] = chr_locus
#print(dicg)
file = open(goin,"r")
fileline = file.readlines()
file.close()
print("Chr\tPos\tEffect\tPvalue\tGenes")
for i in fileline:
    i = i.split()
    chrt = i[0]
    chrtst = int(i[1])-ranges 
    chrted = int(i[1])+ranges 
    headstr = ""
    #print(chrt,chrtst,chrted)
    for k in i:
        headstr += k+"\t"
    headstr += "gene"    
    print(headstr,end = ":")
    for j in  dicg.keys():

        #print(dicg[j])
        gst = int(dicg[j].split("_")[1])
        ged = int(dicg[j].split("_")[2])
        chrg = dicg[j].split("_")[0]
        #print(chrg,gst,ged)
        #break        
        if chrg == chrt and gst > chrtst and ged<chrted:
            print(j,end = "\t")
    print()


基于3V結(jié)果的一個(gè)峰拆分輸出獨(dú)特的QTNs
xlsx轉(zhuǎn)換為csv
xlsx2csv.py

import sys 
import pandas as pd

goin = sys.argv[1]
goout = sys.argv[2]
tblnum = int(sys.argv[3])
def xlsx_to_csv_pd():
    data_xls = pd.read_excel(goin, sheet_name=tblnum)
    data_xls.to_csv(goout, encoding='utf-8')

if __name__ == '__main__':
    xlsx_to_csv_pd()

獨(dú)特QTNs輸出
getposQTNs3V.py

import os
import sys
from numpy import median
goin =  sys.argv[1]
ranges = int(sys.argv[2])
file = open(goin,"r")
lines = list(file.readlines())
dic = {}
dicout = {}
for i in lines:
    i = i.strip()
    if i.find("ID") == -1:
        i = i.replace('"','')
        i = i.split(",")
        chrn = i[4] 
        #print(i)
        dicout[i[4]+"-"+str(i[5])] = [i[6],i[11]]
        if chrn not in dic.keys():
            dic[chrn] = []
            dic[chrn].append(int(i[5]))
        else:
            dic[chrn].append(int(i[5]))



#print(dic)
dicoutmid = {}            
for i in dic.keys():
    
    dicoutmid[i] =[]
    listpos = sorted(dic[i])#,reverse=True)
    #print(listpos)
    for j in range(len(listpos)):
        if j == 0:
            gl = []
            gl.append(listpos[j])

        if j + 1 < len(listpos):
            if listpos[j+1] - listpos[j] < ranges:
                gl.append(listpos[j+1])
            else:
                dicoutmid[i].append(gl)
                gl = []
                gl.append(listpos[j+1])
        else:
            dicoutmid[i].append(gl)
#print(dicoutmid)                
for i in dicoutmid.keys():
    for j in dicoutmid[i]:
        if len(j) %2 == 0:
            j = j[:-1]
        pos = median(j)
        effect = dicout[i+"-"+str(int(round(pos)))][0]
        pvalue = dicout[i+"-"+str(int(round(pos)))][1]
        print(i,int(round(pos)),effect,pvalue)

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容