R語言實現對基因組SNV進行注釋

2022-06-10 03:12:13 字數 1755 閱讀 5822

很多時候,我們需要對取出的snv進行注釋,這個時候可能會在r上進行注釋,通常注釋檔案都含有chr(染色體)、start(開始位點)、end(結束位點)、description(描述),而我們的snv檔案通常是擁有position(位置),因此我們可以先定位chr,再用postion去定位到start和end之間,找到相對應的description。為了加快速度,可以使用二分查詢法。

1

for (value in

dt$value) else

if (value < df[mid,1

])

15if(high18 mid=(low+high)%/%219}

20mid21}

22 }

在r中使用for迴圈效率低,因此也可以用data.table包的foverlap函式,改進**如下,對bed檔案進行注釋,如果要對snv進行注釋,只需要將snv改成相應的start和end相等的bed檔案即可。

1 #!/bin/rscript23

library(data.table)

45 arg <-commandargs(t)

6if (length(arg) != 3

) 14

15 bedfile <- arg[1

]16 annofile <- arg[2

]17 outfile <- arg[3]18

19#read file

20 anno <- fread(annofile,sep="

\t",header=f)

21 bed <- fread(bedfile,sep="

\t",header=f)

22 setnames(anno,c("

v1","

v2","

v3","

v4","

v5","

v9"),c("

chr","

gene

","type

","start

","end

","info"))

23 anno <- anno[type=="

gene

;")[3][[1]],"

\""),function(x)x[2

]))]

24setkey(anno,chr,start,end)

25setkey(bed,v1,v2,v3)

2627

#find overlaps by chr

28 lst <-list()

29for (chri in

intersect(unique(bed$v1),unique(anno$chr)))38}

39 merge_dt <-rbindlist(lst)

40 setnames(merge_dt,c("

v2","

v3","

v4"),c("

start

","end

","name"))

4142

#if one region has more than one gene

43 torm <-list()

44for (i in

1:(nrow(merge_dt)-1))}

45 torm <-unlist(torm)

46 merge_dt <- merge_dt[-torm,]

4748 fwrite(merge_dt,file=outfile)

談談我對基因組區塊鏈的看法

昨天,生信技能樹推送了一篇關於基因區塊鏈的文章。剛好,我對基因測序產品和區塊鏈都很感興趣,所以說幾點我的看法。如果不從區塊鏈屬性來看,僅僅認為它在賣測序服務,4999對於90g的全基因組測序服務而言,也是非常划算。目前公司給我們的 都是65元 g,而dna建庫也要200 300,如果你自己提dna,...

R語言實現KNN 演算法

knn是機器學習中最簡單的分類演算法之一 就是把每乙個測試樣本跟訓練樣本中的每乙個樣本求他們的歐式距離,然後選出最小的幾個,裡面哪乙個類多 這個測試樣本就屬於哪乙個類 用r語言自帶的iris 寫了一下 data iris length iris 1 idx sample 150,100 train ...

R語言實現RMF模型

rmf模型說明 rmf模型是客戶管理中,常被用來衡量客戶價值和客戶創利能力的重要方法。它主要考量三個指標 最近一次消費 recency 近期購買的客戶傾向於再度購買 消費頻率 frequency 經常購買的客戶再次購買概率高 消費金額 monetary 消費金額較多的客戶再次消費可能性更大 根據上述...