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
|
## ----code, echo = FALSE-------------------------------------------------------
code <- function(...) {
cat(paste(..., sep = "\n"))
}
date = "`r doc_date()`"
pkg = "`r pkg_ver('BiocStyle')`"
## ----lastz, eval=FALSE, echo=TRUE---------------------------------------------
# ## lastz aligner
# assemblyDir <- "/Users/gtan/OneDrive/Project/CSC/CNEr/2bit"
# axtDir <- "/Users/gtan/OneDrive/Project/CSC/CNEr/axt"
# assemblyTarget <- file.path(system.file("extdata",
# package="BSgenome.Drerio.UCSC.danRer10"),
# "single_sequences.2bit")
# assemblyQuery <- file.path(system.file("extdata",
# package="BSgenome.Hsapiens.UCSC.hg38"),
# "single_sequences.2bit")
# lavs <- lastz(assemblyTarget, assemblyQuery,
# outputDir=axtDir,
# chrsTarget=c("chr1", "chr2", "chr3"),
# chrsQuery=c("chr1", "chr2", "chr3"),
# distance="far", mc.cores=4)
#
# ## lav files to psl files conversion
# psls <- lavToPsl(lavs, removeLav=FALSE, binary="lavToPsl")
## ----last, eval=FALSE, echo=TRUE----------------------------------------------
# ## Build the lastdb index
# system2(command="lastdb", args=c("-c", file.path(assemblyDir, "danRer10"),
# file.path(assemblyDir, "danRer10.fa")))
#
# ## Run last aligner
# lastal(db=file.path(assemblyDir, "danRer10"),
# queryFn=file.path(assemblyDir, "hg38.fa"),
# outputFn=file.path(axtDir, "danRer10.hg38.maf"),
# distance="far", binary="lastal", mc.cores=4L)
#
# ## maf to psl
# psls <- file.path(axtDir, "danRer10.hg38.psl")
# system2(command="maf-convert", args=c("psl",
# file.path(axtDir, "danRer10.hg38.maf"),
# ">", psls))
## ----chain, eval=FALSE, echo=TRUE---------------------------------------------
# ## Join close alignments
# chains <- axtChain(psls, assemblyTarget=assemblyTarget,
# assemblyQuery=assemblyQuery, distance="far",
# removePsl=FALSE, binary="axtChain")
#
# ## Sort and combine
# allChain <- chainMergeSort(chains, assemblyTarget, assemblyQuery,
# allChain=file.path(axtDir,
# paste0(sub("\\.2bit$", "", basename(assemblyTarget),
# ignore.case=TRUE), ".",
# sub("\\.2bit$", "", basename(assemblyQuery),
# ignore.case=TRUE), ".all.chain")),
# removeChains=FALSE, binary="chainMergeSort")
## ----netting, eval=FALSE, echo=TRUE-------------------------------------------
# ## Filtering out chains
# allPreChain <- chainPreNet(allChain, assemblyTarget, assemblyQuery,
# allPreChain=file.path(axtDir,
# paste0(sub("\\.2bit$", "",
# basename(assemblyTarget),
# ignore.case = TRUE), ".",
# sub("\\.2bit$", "",
# basename(assemblyQuery),
# ignore.case = TRUE),
# ".all.pre.chain")),
# removeAllChain=FALSE, binary="chainPreNet")
#
# ## Keep the best chain and add synteny information
# netSyntenicFile <- chainNetSyntenic(allPreChain, assemblyTarget, assemblyQuery,
# netSyntenicFile=file.path(axtDir,
# paste0(sub("\\.2bit$", "",
# basename(assemblyTarget),
# ignore.case = TRUE), ".",
# sub("\\.2bit$", "",
# basename(assemblyQuery),
# ignore.case = TRUE),
# ".noClass.net")),
# binaryChainNet="chainNet", binaryNetSyntenic="netSyntenic")
## ----axtNet, eval=FALSE, echo=TRUE--------------------------------------------
# netToAxt(netSyntenicFile, allPreChain, assemblyTarget, assemblyQuery,
# axtFile=file.path(axtDir,
# paste0(sub("\\.2bit$", "",
# basename(assemblyTarget),
# ignore.case = TRUE), ".",
# sub("\\.2bit$", "",
# basename(assemblyQuery),
# ignore.case = TRUE),
# ".net.axt")),
# removeFiles=FALSE,
# binaryNetToAxt="netToAxt", binaryAxtSort="axtSort")
#
## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()
|
