Pipeline of generating pairwise whole genome alignment.
CNEr 1.14.0
CNEr identifies conserved noncoding elments (CNEs) from pairwise whole genome alignment (net.axt files) of two species. UCSC has provided alignments between many species on the downloads, hence it is highly recommended to use their alignments when available. When the alignments of some new assemblies/species are not availble from UCSC yet, this vignette describes the pipeline of generating the alignments merely from soft-masked 2bit files or fasta files. This vignette is based on the genomewiki from UCSC.
Note:
List of external softwares have to be installed on the machine: * The sequence alignment program LASTZ * Kent Utilities. In this pipeline, lavToPsl, axtChain, chainMergeSort, chainPreNet, chainNet, netSyntenic, netToAxt, axtSort are essential. netClass is optional.
Here, as an example, we will get the pairwise alignment only on “chr1”, “chr2” and “chr3” between zebrafish(danRer10) and human(hg38).
First we need to download the 2bit files from UCSC, and set the appropriate paths of assemblyTarget, assemblyQuery and intermediate files. Then we can run lastz function to generate the lav files.
## 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")One essential argument here is the distance. It determines the scoring matrix to use in lastz aligner. See ?scoringMatrix for more details.
Note: This step may encounter difficulties if the two assemblies are too fragmented, because there can be millions of combinations of the chromosomes/scaffolds. The lastz is overwhelmed with reading small pieces from assembly for each combination, rather than doing actual alignment. In this case, another aligner last is recommended and introduced in nex section.
last aligner is considered faster and memory efficient. It creates maf file, which can converted to psl files. Then the same following processes can be used on psl files.
Different from lastz, last aligner starts with fasta files. The target genome sequence has to build the index file first, and then align with the query genome sequence.
## 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))Another alternative of alignment software is YASS. It may be added into this pipeline after we test the performance.
If two matching alignments next to each other are close enough, they are joined into one fragment. Then these chain files are sorted and combined into one big file.
## 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")In this step, first we filter out chains that are unlikely to be netted by chainPreNet. During the alignment, every genomic fragment can match with several others, and certainly we want to keep the best one. This is done by chainNet. Then we add the synteny information with netSyntenic.
## 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")As the last step, we create the .net.axt file from the previous net and chain files.
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")Here is the output of sessionInfo() on the system on which this document was compiled:
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] Gviz_1.22.0 BSgenome.Ggallus.UCSC.galGal3_1.4.0
## [3] BSgenome.Hsapiens.UCSC.hg19_1.4.0 BSgenome_1.46.0
## [5] rtracklayer_1.38.0 Biostrings_2.46.0
## [7] XVector_0.18.0 GenomicRanges_1.30.0
## [9] GenomeInfoDb_1.14.0 IRanges_2.12.0
## [11] S4Vectors_0.16.0 BiocGenerics_0.24.0
## [13] CNEr_1.14.0 BiocStyle_2.6.0
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.10.0 bitops_1.0-6
## [3] matrixStats_0.52.2 bit64_0.9-7
## [5] RColorBrewer_1.1-2 progress_1.1.2
## [7] httr_1.3.1 rprojroot_1.2
## [9] tools_3.4.2 backports_1.1.1
## [11] R6_2.2.2 rpart_4.1-11
## [13] Hmisc_4.0-3 DBI_0.7
## [15] lazyeval_0.2.1 colorspace_1.3-2
## [17] nnet_7.3-12 RMySQL_0.10.13
## [19] gridExtra_2.3 prettyunits_1.0.2
## [21] curl_3.0 bit_1.1-12
## [23] compiler_3.4.2 Biobase_2.38.0
## [25] htmlTable_1.9 DelayedArray_0.4.0
## [27] labeling_0.3 bookdown_0.5
## [29] checkmate_1.8.5 scales_0.5.0
## [31] readr_1.1.1 stringr_1.2.0
## [33] digest_0.6.12 Rsamtools_1.30.0
## [35] foreign_0.8-69 rmarkdown_1.6
## [37] R.utils_2.5.0 dichromat_2.0-0
## [39] base64enc_0.1-3 pkgconfig_2.0.1
## [41] htmltools_0.3.6 ensembldb_2.2.0
## [43] htmlwidgets_0.9 rlang_0.1.2
## [45] RSQLite_2.0 VGAM_1.0-4
## [47] BiocInstaller_1.28.0 shiny_1.0.5
## [49] BiocParallel_1.12.0 acepack_1.4.1
## [51] R.oo_1.21.0 VariantAnnotation_1.24.0
## [53] RCurl_1.95-4.8 magrittr_1.5
## [55] GO.db_3.4.2 GenomeInfoDbData_0.99.1
## [57] Formula_1.2-2 Matrix_1.2-11
## [59] Rcpp_0.12.13 munsell_0.4.3
## [61] R.methodsS3_1.7.1 stringi_1.1.5
## [63] yaml_2.1.14 SummarizedExperiment_1.8.0
## [65] zlibbioc_1.24.0 AnnotationHub_2.10.0
## [67] plyr_1.8.4 blob_1.1.0
## [69] lattice_0.20-35 splines_3.4.2
## [71] GenomicFeatures_1.30.0 annotate_1.56.0
## [73] hms_0.3 KEGGREST_1.18.0
## [75] knitr_1.17 reshape2_1.4.2
## [77] biomaRt_2.34.0 XML_3.98-1.9
## [79] evaluate_0.10.1 biovizBase_1.26.0
## [81] latticeExtra_0.6-28 data.table_1.10.4-3
## [83] httpuv_1.3.5 png_0.1-7
## [85] gtable_0.2.0 poweRlaw_0.70.1
## [87] assertthat_0.2.0 ggplot2_2.2.1
## [89] mime_0.5 xtable_1.8-2
## [91] AnnotationFilter_1.2.0 survival_2.41-3
## [93] tibble_1.3.4 GenomicAlignments_1.14.0
## [95] AnnotationDbi_1.40.0 memoise_1.1.0
## [97] cluster_2.0.6 interactiveDisplayBase_1.16.0