-
Notifications
You must be signed in to change notification settings - Fork 0
/
module4_visualization.Rmd
executable file
·74 lines (66 loc) · 2.83 KB
/
module4_visualization.Rmd
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
---
title: "module4_visualization"
date: "1/31/2022"
output: html_document
---
Load data
```{r}
ratio_dataTable <- read.table(file = "https://github.com/bioinformatics-ca/CAN_2021/raw/main/Module4/Data/CBW_regions_c0098_Tumor.sorted.markduplicates.bam_sample.cpn_ratio.txt",header=TRUE)
ratio <- data.frame(ratio_dataTable)
BAF_dataTable <- read.table(file = "https://github.com/bioinformatics-ca/CAN_2021/raw/main/Module4/Data/CBW_regions_c0098_Tumor.sorted.markduplicates.bam_sample.cpn_BAF.txt", header=TRUE)
BAF<-data.frame(BAF_dataTable)
ploidy <- 2
```
Plot
```{r}
for (chrom in c(3))
{
tt <- which(ratio$Chromosome==chrom)
if (length(tt)>0)
{
plot(ratio$Start[tt],log2(ratio$Ratio[tt]),xlab = paste ("position, chr",chrom),ylab = "normalized copy number profile (log2)",pch = ".",col = colors()[88],cex=4)
tt <- which(ratio$Chromosome==chrom & ratio$CopyNumber>ploidy )
points(ratio$Start[tt],log2(ratio$Ratio[tt]),pch = ".",col = colors()[136],cex=4)
tt <- which(ratio$Chromosome==chrom & ratio$CopyNumber<ploidy & ratio$CopyNumber!= -1)
points(ratio$Start[tt],log2(ratio$Ratio[tt]),pch = ".",col = colors()[461],cex=4)
tt <- which(ratio$Chromosome==chrom)
}
tt <- which(ratio$Chromosome==chrom)
}
```
BAF plot
```{r}
for (i in c(3))
{
tt <- which(BAF$Chromosome==i)
if (length(tt)>0)
{
lBAF <-BAF[tt,]
plot(lBAF$Position,lBAF$BAF,ylim = c(-0.1,1.1),xlab = paste ("position, chr",i),ylab = "BAF",pch = ".",col = colors()[1])
tt <- which(lBAF$A==0.5)
points(lBAF$Position[tt],lBAF$BAF[tt],pch = ".",col = colors()[92])
tt <- which(lBAF$A!=0.5 & lBAF$A>=0)
points(lBAF$Position[tt],lBAF$BAF[tt],pch = ".",col = colors()[62])
tt <- 1
pres <- 1
if (length(lBAF$A)>4) {
for (j in c(2:(length(lBAF$A)-pres-1))) {
if (lBAF$A[j]==lBAF$A[j+pres]) {
tt[length(tt)+1] <- j
}
}
points(lBAF$Position[tt],lBAF$A[tt],pch = ".",col = colors()[24],cex=3)
points(lBAF$Position[tt],lBAF$B[tt],pch = ".",col = colors()[24],cex=3)
}
tt <- 1
pres <- 1
if (length(lBAF$FittedA)>4) {
for (j in c(2:(length(lBAF$FittedA)-pres-1))) {
if (lBAF$FittedA[j]==lBAF$FittedA[j+pres]) {
tt[length(tt)+1] <- j
}
}
points(lBAF$Position[tt],lBAF$FittedA[tt],pch = ".",col = colors()[463],cex=3)
points(lBAF$Position[tt],lBAF$FittedB[tt],pch = ".",col = colors()[463],cex=3)
}}}
```