-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscript_use_StruturalVariantAnnotations.R
108 lines (80 loc) · 3.85 KB
/
script_use_StruturalVariantAnnotations.R
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
##########################################################
##########################################################
# CRAN packages
library(devtools)
library(stringr)
library(dplyr)
library(ggplot2)
# bioconductor packages
library(VariantAnnotation)
library(GenomicRanges)
library(GenomicFeatures)
library(rtracklayer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
##########################################################
##########################################################
library(StructuralVariantAnnotation)
# install_github("d-cameron/StructuralVariantAnnotation")
##########################################################
##########################################################
# a <- "./vcf.DELLY.SPCG-EWS008_2D.vcf"
a <- "./vcf.DELLY.SPCG-EWS008_13R.vcf"
vcf <- readVcf(a, "hg38")
rowRanges(vcf)
info(vcf)
info(header(vcf))
geno(vcf)
geno(header(vcf))
##########################################################
# breakpointRanges
# extracts the structural variants as a BreakendGRanges
##########################################################
# convert to breakend GRanges
##########################################################
##########################################################
gr <- breakpointRanges(vcf)
# retain only primary chromosomes
seqlevelsStyle(gr) <- "UCSC"
gr <- gr[seqnames(gr) %in% paste0("chr", c(1:22, "X", "Y"))]
seqlevels(gr) <- paste0("chr", c(1:22, "X", "Y"))
gr <- breakpointRanges(vcf)
##########################################################
##########################################################
gr
# remove breakends that now don't have a partner (eg: chr1 -> chrMT)
gr <- gr[gr$partner %in% names(gr)]
# annotate breakends with gene names and gene orientation
gns <- genes(TxDb.Hsapiens.UCSC.hg38.knownGene)
hits <- as.data.frame(findOverlaps(gr, gns, ignore.strand=TRUE))
hits$SYMBOL <- biomaRt::select(org.Hs.eg.db, gns[hits$subjectHits]$gene_id, "SYMBOL")$SYMBOL
hits$gene_strand <- as.character(strand(gns[hits$subjectHits]))
hits <- hits %>%
group_by(queryHits) %>%
summarise(SYMBOL=paste(SYMBOL, collapse=","), gene_strand=paste0(gene_strand, collapse=""))
gr$SYMBOL <- ""
gr$geneStrand <- ""
gr$SYMBOL[hits$queryHits] <- hits$SYMBOL
gr$geneStrand[hits$queryHits] <- hits$gene_strand
# require the breakpoint to be between different genes
gr <- gr[gr$SYMBOL != partner(gr)$SYMBOL,]
gr <- gr[gr$SYMBOL != "" & partner(gr)$SYMBOL != "",]
##########################################################
##########################################################
##########################################################
# requires the breakpoint to possibly generate a fusion transcript
##########################################################
##########################################################
##########################################################
gr$couldBeThreePrimeStart <- str_detect(gr$geneStrand, stringr::fixed(as.character(strand(gr))))
gr$couldBeFivePrimeEnd <- str_detect(gr$geneStrand, stringr::fixed(ifelse(as.character(strand(gr))=="+", "-", "+")))
gr <- gr[(gr$couldBeThreePrimeStart & partner(gr)$couldBeFivePrimeEnd) |
(gr$couldBeFivePrimeEnd & partner(gr)$couldBeThreePrimeStart),]
##########################################################
##########################################################
########################################################## writing the output :
write.table(gr, file="output_StructuralVariantAnnotation",
sep="\t", row.names=F)
####################################################################################################################################
####################################################################################################################################
####################################################################################################################################