scripts/cn.mops.R at master · PBR/scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
dir()
getwd()
setwd('/media/bulk_01/users/goryu003/samples')
dir()
ls()
dir()
df_1_size=read.table("1_size_minwidth_cnv.txt", as.is=T)
bed=read.table("ITAG2.3_gene_models_ch04_0-5M-3.bed", as.is=T)
head(df_1_size)
df_1_size[1,3]
head(bed)
bed[1,2]
dir()
library(cn.mops)
BAMFiles <- list.files(pattern=".bam$")
str(BAMFiles)
is.list(BAMFiles)
is.vector(BAMFiles)
BAMFiles
names=c("Moneymaker","S. Arcanum","S. habrochaites f. glabratum","S. pennellii")
segments <- read.table("ITAG2.3_gene_models_ch04_0-5M-3.bed",sep="\t",as.is=TRUE)
segments
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
gr
ls()
dir()
count <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,mode="paired",sampleNames=names)
resCNMOPS <- cn.mops(count, minWidth=1)
resCNMOPS <- calcIntegerCopyNumbers(resCNMOPS)
resCNMOPS
help(CNVDetectionResult)
integerCopyNumber(resCNMOPS)
as.data.frame(integerCopyNumber(resCNMOPS))
count_table=as.data.frame(integerCopyNumber(resCNMOPS))
head(count_table)
head(count_table)
write.table(count_table, "count_table.txt", row.names=F)
head(count_table[,c(1:4,6:9])
head(count_table[,c(1:4,6:9)])
count_table_2=count_table[,c(1:4,6:9)]
count_table_2
head(count_table_2)
write.table(count_table_2, "count_table.txt", row.names=F)
dir()
source("subset_CNV_predict.R")
head(count_table_2)
head(count_table_2[,5])
head(count_table_2[,5:8])
source("subset_CNV_predict.R")
ls()
bed
ls()
dir()
CNV=read.table("named_cnv.txt",as.is=T)
CNV_final = merge(CNV,bed,by.x=c(1,2,3),by.y=c(1,2,3),all.x=T,all.y=F)
head(CNV_final)
write.table(CNV_final,"CNV_final.txt",row.names=F)
write.table(CNV_final,"CNV_final.txt",row.names=F,col.names=F)