Perform alternative splicing analysis for multiple Visium samples

This vignette demonstrates how to analyze alternative splicing across multiple Visium samples. The procedure for spatial transcriptomics follows the same principles as scRNA-seq but requires careful handling of imaging data and spatial coordinates for accurate visualization. Additionally, this example illustrates how to visualize genome tracks from multiple BAM files.

0. Prepare the data.

We will use four Visium slices for 10X data source page. These slices come from two sections, each section contain two slices, one is sagittal anterior, and another one is sagittal posterior of mouse brain. We download the processed BAM files, and annotated them with PISA.

Section 1, sagittal anterior

Section 1, sagittal posterior

Section 2, sagittal anterior

Section 2, sagittal posterior

The following commands run on a terminal.

# Create directory for each section and slice
mkdir -p visium/section1/Posterior
mkdir -p visium/section1/Anterior
mkdir -p visium/section2/Posterior
mkdir -p visium/section2/Anterior

# Enter to work directory and download GTF for mouse
cd visium;
wget -c https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M10/gencode.vM10.annotation.gtf.gz

# Download the bam file and image files to each directory
cd section1/Posterior
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_possorted_genome_bam.bam
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_possorted_genome_bam.bam.bai
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_spatial.tar.gz

# unpack image files
tar zxvf V1_Mouse_Brain_Sagittal_Posterior_spatial.tar.gz

# Annotate and count features
PISA anno -gtf ../../gencode.vM10.annotation.gtf.gz -exon -psi -o anno.bam -t 7 V1_Mouse_Brain_Sagittal_Posterior_possorted_genome_bam.bam

mkdir exp
mkdir exon
mkdir junction
mkdir exclude
PISA count -tag CB -umi UB -anno-tag GN -outdir exp anno.bam
PISA count -tag CB -umi UB -anno-tag EX -outdir exon anno.bam
PISA count -tag CB -umi UB -anno-tag JC -outdir junction anno.bam
PISA count -tag CB -umi UB -anno-tag ER -outdir exclude anno.bam

## Go back to work directory, and reenter to another slice's directory
cd -; 
cd section1/Anterior
wget -c https://s3-us-west-2.amazonaws.com/10x.files/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_possorted_genome_bam.bam
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_possorted_genome_bam.bam.bai
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz

# unpack image file
tar zxvf V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz

# Annotate and count features
PISA anno -gtf ../../gencode.vM10.annotation.gtf.gz -exon -psi -o anno.bam -t 7 V1_Mouse_Brain_Sagittal_Anterior_possorted_genome_bam.bam

mkdir exp
mkdir exon
mkdir junction
mkdir exclude
PISA count -tag CB -umi UB -anno-tag GN -outdir exp anno.bam
PISA count -tag CB -umi UB -anno-tag EX -outdir exon anno.bam
PISA count -tag CB -umi UB -anno-tag JC -outdir junction anno.bam
PISA count -tag CB -umi UB -anno-tag ER -outdir exclude anno.bam

# Swith to section 2 
cd -;
cd section2/Anterior
wget -c https://s3-us-west-2.amazonaws.com/10x.files/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior_Section_2/V1_Mouse_Brain_Sagittal_Anterior_Section_2_possorted_genome_bam.bam
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior_Section_2/V1_Mouse_Brain_Sagittal_Anterior_Section_2_possorted_genome_bam.bam.bai
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior_Section_2/V1_Mouse_Brain_Sagittal_Anterior_Section_2_spatial.tar.gz

tar zxvf V1_Mouse_Brain_Sagittal_Anterior_Section_2_spatial.tar.gz

PISA anno -gtf ../../gencode.vM10.annotation.gtf.gz -exon -psi -o anno.bam -t 7 V1_Mouse_Brain_Sagittal_Anterior_Section_2_possorted_genome_bam.bam

mkdir exp
mkdir exon
mkdir junction
mkdir exclude
PISA count -tag CB -umi UB -anno-tag GN -outdir exp anno.bam
PISA count -tag CB -umi UB -anno-tag EX -outdir exon anno.bam
PISA count -tag CB -umi UB -anno-tag JC -outdir junction anno.bam
PISA count -tag CB -umi UB -anno-tag ER -outdir exclude anno.bam


cd -;
cd section2/Sagittal
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior_Section_2/V1_Mouse_Brain_Sagittal_Posterior_Section_2_possorted_genome_bam.bam
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior_Section_2/V1_Mouse_Brain_Sagittal_Posterior_Section_2_possorted_genome_bam.bam.bai
wget -c https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior_Section_2/V1_Mouse_Brain_Sagittal_Posterior_Section_2_spatial.tar.gz

tar zxvf V1_Mouse_Brain_Sagittal_Posterior_Section_2_spatial.tar.gz

PISA anno -gtf ../../gencode.vM10.annotation.gtf.gz -exon -psi -o anno.bam -t 7 V1_Mouse_Brain_Sagittal_Posterior_Section_2_possorted_genome_bam.bam

mkdir exp
mkdir exon
mkdir junction
mkdir exclude
PISA count -tag CB -umi UB -anno-tag GN -outdir exp anno.bam
PISA count -tag CB -umi UB -anno-tag EX -outdir exon anno.bam
PISA count -tag CB -umi UB -anno-tag JC -outdir junction anno.bam
PISA count -tag CB -umi UB -anno-tag ER -outdir exclude anno.bam

# Now go back to work directory
cd ../../..

1. Preprocess all four slices and intergrate them into one object

This section will perform Seurat pipeline for all the sections. A Seurat workflow for one section can be found here https://satijalab.org/seurat/articles/spatial_vignette. The difference between this workflow with Seurat’s is we use counts generated from PISA, while Seurat use counts generated with CellRanger’s pipeline. For the details of different between PISA and CellRanger, please see our manual.

require(Yano)
Loading required package: Yano
── Attaching packages ────────────────────────────────────────────── Yano 1.0 ──
✔ dplyr   1.1.4     ✔ Seurat  5.1.0
✔ ggplot2 3.5.1     
exp.1 <- ReadPISA("./visium/section1/Anterior/exp/")
image.1 <- Read10X_Image("./visium/section1/Anterior/spatial/")

## QuickRecipe() is actually an integration function of Seurat workflow
obj.1 <- QuickRecipe(exp.1[, Cells(image.1)], verbose = FALSE)
Normalizing layer: counts
obj.1[['slice1']] <- image.1[colnames(obj.1),]
obj.1$orig.ident <- "sec1_ant"


exp.2 <- ReadPISA("./visium/section1/Posterior/exp/")
image.2 <- Read10X_Image("./visium/section1/Posterior/spatial/")
obj.2 <- QuickRecipe(exp.2[,Cells(image.2)], verbose = FALSE)
Normalizing layer: counts
obj.2[['slice2']] <- image.2[colnames(obj.2),]
obj.2$orig.ident <- "sec1_pos"

exp.3 <- ReadPISA("./visium/section2/Anterior/exp/")
image.3 <- Read10X_Image("./visium/section2/Anterior/spatial/")
obj.3 <- QuickRecipe(exp.3[, Cells(image.3)], verbose = FALSE)
Normalizing layer: counts
obj.3[['slice3']] <- image.3[colnames(obj.3),]
obj.3$orig.ident <- "sec2_ant"

exp.4 <- ReadPISA("./visium/section2/Posterior/exp/")
image.4 <- Read10X_Image("./visium/section2/Posterior/spatial/")
obj.4 <- QuickRecipe(exp.4[, Cells(image.4)], verbose = FALSE)
Normalizing layer: counts
obj.4[['slice4']] <- image.4[colnames(obj.4),]
obj.4$orig.ident <- "sec2_pos"

# Merge all processed Seurat object
obj <- merge(obj.1, y = list(obj.2, obj.3, obj.4), add.cell.ids = c("sec1_ant", "sec1_pos", "sec2_ant", "sec2_pos"))

obj <- QuickRecipe(obj)
Set default assay to RNA
object <- NormalizeData(object, normalization.method = "LogNormalize", scale.factor = 10000)
Normalizing layer: counts.1
Normalizing layer: counts.2
Normalizing layer: counts.3
Normalizing layer: counts.4
object <- FindVariableFeatures(object, selection.method = "vst", nfeatures = 3000)
Finding variable features for layer counts.1
Finding variable features for layer counts.2
Finding variable features for layer counts.3
Finding variable features for layer counts.4
object <- ScaleData(object, features =  @features, vars.to.regress = "nCount_RNA")
Regressing out nCount_RNA
Centering and scaling data matrix
object <- RunPCA(object, features = @features
PC_ 1 
Positive:  Dbi, Pcp2, Car8, Sparc, Zic1, Cbln1, Pvalb, Ppp1r17, Lhx1os, Homer3 
       Sept4, Metrn, Gm2694, Atp2a3, Hopx, Inpp5a, Ank1, Fam107b, Neurod1, Kcnc3 
       Timp4, Icmt, Col18a1, Aldoc, Tspan9, Dner, Slc1a6, Calb2, Chn2, Kcng4 
Negative:  Nrgn, Ctxn1, Cx3cl1, Chn1, Ngef, Ptk2b, Camkv, Camk2n1, Pde2a, Hpca 
       Arpp21, Fbxl16, Enc1, Ddn, Pld3, Basp1, Ppp3ca, Rgs4, Chst1, Camk2a 
       Arf3, Chrm1, Ptpn5, Icam5, Phactr1, Hpcal4, Rprml, Kcnip2, Syt5, Foxg1 
PC_ 2 
Positive:  Snap25, Neurod1, Snhg11, Dnm1, Adcy1, Stmn2, Ywhah, Scn1b, Stmn3, Gm2694 
       Cbln1, Cacna1a, Shf, Sh3gl2, Stxbp1, Grin2c, Gm42742, Zbtb18, Nrep, Mir124a-1hg 
       Diras2, Tmem132a, Meg3, Sphkap, Grm4, Ndrg3, Dgkz, Nptx1, Car8, Cadm3 
Negative:  Apod, Plp1, Mal, Cldn11, Gsn, Enpp2, Mbp, Mobp, Cryab, Mog 
       Mag, Tmsb4x, Car2, Gfap, Cnp, Plekhb1, Trf, Ndrg1, Ermn, Ugt8a 
       Efnb3, Bcas1, Gatm, Ppp1r14a, Fa2h, Gpr37, Pllp, Pdlim2, Gltp, Gm20425 
PC_ 3 
Positive:  Penk, Rgs9, Pde10a, Adora2a, Ppp1r1b, Gpr88, Drd1, Tac1, Syndig1l, Gpr6 
       Rxrg, Lrrc10b, Pde1b, Adcy5, Rasd2, Scn4b, Tmem158, Gng7, Rarb, Slc32a1 
       Slc35d3, Actn2, Strip2, Ankrd63, Serpina9, Pcp4l1, Meis2, Sez6, Gnal, Pcp4 
Negative:  Cck, 3110035E14Rik, Nrn1, Olfm1, Stmn1, Slc17a7, Lingo1, Stx1a, Tbr1, Gnas 
       Slc30a3, Vsnl1, Dkk3, Ncald, 1110008P14Rik, Efhd2, Basp1, Rtn4r, Neurod6, Miat 
       Gm11549, Mpped1, Pde1a, Mical2, Fxyd7, Nptxr, E130012A19Rik, Gabbr2, Satb2, Syt13 
PC_ 4 
Positive:  Mbp, Mobp, Qdpr, Plp1, Mag, Trf, Plekhb1, Cnp, Bcas1, Mal 
       Mog, Tspan2, Ermn, Sept4, Cryab, Car2, Cldn11, Nefm, Pllp, Ppp1r14a 
       Tmem88b, Gpr37, Vamp1, Ugt8a, Gatm, Fa2h, Scn4b, Nefl, Nefh, Prr18 
Negative:  Islr, Igf2, Igfbp2, Mgp, Fmod, Slc6a20a, Slc13a4, Ogn, Aebp1, Slc6a13 
       Fabp7, Aldh1a2, Gjb2, Slc22a6, Col1a2, Rbp1, Fn1, Slc13a3, Bmp7, Crabp2 
       Pcolce, Serping1, Efemp1, Col1a1, Ifitm3, Nupr1, Bgn, Cpne6, Nbl1, Thbd 
PC_ 5 
Positive:  Gng4, Slc6a11, Ptpro, Sp8, Dlx1, Nrip3, Hap1, Cdhr1, Scg2, Doc2g 
       Th, Pcbp3, Tubb2a, Cpne4, Dcx, Scgn, Ppm1e, Tshz1, Zcchc12, A230065H16Rik 
       Pbx3, Dlx2, Lgr6, Lgr5, Cacng5, Sp9, Wdr6, Tmem130, Adamts19, Nmb 
Negative:  Itpka, Itpr1, Zbtb18, Hpca, Ppp1r1b, Arpp19, Cplx2, Rnf112, Adcy1, Camk4 
       Cnr1, Nptx1, Mmp17, Cnksr2, Lmo4, Prkcb, Sepw1, Neurod2, Calb1, Lzts3 
       Sptbn2, Tmem132a, Tesc, Foxp1, Fbxl16, Lingo3, Dkk3, Ppp3ca, Atp1a1, Lrrc10b 
object <- FindNeighbors(object, dims = 1:20)
Computing nearest neighbor graph
Computing SNN
object <- FindClusters(object, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 12148
Number of edges: 386668

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9273
Number of communities: 17
Elapsed time: 0 seconds
object <- RunUMAP(object, dims = 1:20)
23:29:59 UMAP embedding parameters a = 0.9922 b = 1.112
23:29:59 Read 12148 rows and found 20 numeric columns
23:29:59 Using Annoy for neighbor search, n_neighbors = 30
23:29:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:30:00 Writing NN index file to temp file /var/folders/38/qfwf72jx3413kwjx7lgqv3xr0000gn/T//RtmpWK6jSu/file85b761505099
23:30:00 Searching Annoy index using 1 thread, search_k = 3000
23:30:02 Annoy recall = 100%
23:30:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
23:30:02 Initializing from normalized Laplacian + noise (using RSpectra)
23:30:02 Commencing optimization for 200 epochs, with 498054 positive edges
23:30:07 Optimization finished
p1 <- DimPlot(obj, group.by = "orig.ident")
p2 <- DimPlot(obj)
p1 + p2

SpatialPlot(obj, pt.size.factor = 2) 

exon.1 <- ReadPISA("./visium/section1/Anterior/exon/", prefix = "sec1_ant_", cells = Cells(obj))
exon.2 <- ReadPISA("./visium/section1/Posterior/exon/", prefix = "sec1_pos_", cells = Cells(obj))
exon.3 <- ReadPISA("./visium/section2/Anterior/exon/", prefix = "sec2_ant_", cells = Cells(obj))
exon.4 <- ReadPISA("./visium/section2/Posterior/exon/", prefix = "sec2_pos_", cells = Cells(obj))
exon <- mergeMatrix(exon.1, exon.2, exon.3, exon.4)
'as(<dgTMatrix>, "dgCMatrix")' is deprecated.
Use 'as(., "CsparseMatrix")' instead.
See help("Deprecated") and help("Matrix-deprecated").
obj[['exon']] <- CreateAssayObject(exon[, Cells(obj)], min.cells = 20)
DefaultAssay(obj) <- "exon"
obj <- NormalizeData(obj)
obj <- ParseExonName(obj)
Working on assay exon
obj <- RunAutoCorr(obj)
Working on assay : exon
Run autocorrelation test for 80178 features.
Runtime : 40.62745 secs
obj <- SetAutoCorrFeatures(obj)
39414 autocorrelated features.
obj <- RunBlockCorr(obj, bind.name = "gene_name", bind.assay = "RNA")
Working on assay exon.
Working on binding assay RNA.
Processing 39414 features.
Processing 15785 blocks.
Retrieve binding data from assay RNA.
Use "data" layer for test features and binding features.
Using 7 threads.
Using method "D" with mode 1.
Use predefined weight matrix "pca_wm".
Runtime : 30.41442 mins.
FbtPlot(obj, val = "gene_name.padj")

obj[['exon']][[]] %>% filter(gene_name.pval<1e-40)
                                   chr     start       end gene_name strand
chr11:30106785-30109152/-/Sptbn1 chr11  30106785  30109152    Sptbn1      -
chr11:53548291-53549565/+/Sept8  chr11  53548291  53549565     Sept8      +
chr12:111806315-111806775/+/Klc1 chr12 111806315 111806775      Klc1      +
chr12:111806315-111806727/+/Klc1 chr12 111806315 111806727      Klc1      +
chr12:111806315-111806711/+/Klc1 chr12 111806315 111806711      Klc1      +
chr12:111806315-111806726/+/Klc1 chr12 111806315 111806726      Klc1      +
chr17:25717528-25717581/+/Gng13  chr17  25717528  25717581     Gng13      +
chrX:163926903-163929546/+/Ap1s2  chrX 163926903 163929546     Ap1s2      +
                                 moransi.pval    moransi autocorr.variable
chr11:30106785-30109152/-/Sptbn1            0 0.15018334              TRUE
chr11:53548291-53549565/+/Sept8             0 0.41877021              TRUE
chr12:111806315-111806775/+/Klc1            0 0.10144477              TRUE
chr12:111806315-111806727/+/Klc1            0 0.10131377              TRUE
chr12:111806315-111806711/+/Klc1            0 0.09578263              TRUE
chr12:111806315-111806726/+/Klc1            0 0.09829769              TRUE
chr17:25717528-25717581/+/Gng13             0 0.15957157              TRUE
chrX:163926903-163929546/+/Ap1s2            0 0.14138305              TRUE
                                 gene_name.D gene_name.r gene_name.pval
chr11:30106785-30109152/-/Sptbn1   0.4957796 -0.02554239   7.057288e-60
chr11:53548291-53549565/+/Sept8    0.3818213  0.05181663   1.718951e-44
chr12:111806315-111806775/+/Klc1   0.4438596  0.09998905   2.876458e-55
chr12:111806315-111806727/+/Klc1   0.4439723  0.10056806   2.817245e-55
chr12:111806315-111806711/+/Klc1   0.4384527  0.11168103   1.359610e-53
chr12:111806315-111806726/+/Klc1   0.4388710  0.10047243   1.271424e-53
chr17:25717528-25717581/+/Gng13    0.3935489 -0.02563280   3.712364e-43
chrX:163926903-163929546/+/Ap1s2   0.3609878  0.02473633   4.911901e-46
                                 gene_name.mean gene_name.var gene_name.padj
chr11:30106785-30109152/-/Sptbn1      0.1250099   0.010071091   2.735405e-55
chr11:53548291-53549565/+/Sept8       0.1223685   0.010483515   9.518079e-41
chr12:111806315-111806775/+/Klc1      0.1257658   0.009702447   3.716384e-51
chr12:111806315-111806727/+/Klc1      0.1257633   0.009703739   3.716384e-51
chr12:111806315-111806711/+/Klc1      0.1254133   0.009964339   1.053970e-49
chr12:111806315-111806726/+/Klc1      0.1258574   0.009956102   1.053970e-49
chr17:25717528-25717581/+/Gng13       0.1232859   0.011322626   1.798640e-39
chrX:163926903-163929546/+/Ap1s2      0.1235875   0.009203270   3.173088e-42
SpatialPlot(obj, features = c("chr11:53548291-53549565/+/Sept8","Sept8"), pt.size.factor = 2)

2. Visualize genome tracks of multiple samples

Since the spot/cell names have already been renamed in the Seurat object, we need to define a binding list that maps the new cell names to their original names, which will be used for reading sequences from the BAM files. In this vignette, we use four BAM files, requiring us to define four corresponding cell vectors. Each vector contains group factor values, where the names of the values correspond to the original cell names in the BAM file.

gtf <- gtf2db("./visium/gencode.vM10.annotation.gtf.gz")
[2025-02-20 00:02:14] GTF loading..
[2025-02-20 00:02:20] Load 48440 genes.
bamfiles <- list(
  sec1_ant = "./visium/section1/Anterior/V1_Mouse_Brain_Sagittal_Anterior_possorted_genome_bam.bam",
  sec1_pos = "./visium/section1/Posterior/V1_Mouse_Brain_Sagittal_Posterior_possorted_genome_bam.bam",
  sec2_ant = "./visium/section2/Anterior/V1_Mouse_Brain_Sagittal_Anterior_Section_2_possorted_genome_bam.bam",
  sec2_pos = "./visium/section2/Posterior/V1_Mouse_Brain_Sagittal_Posterior_Section_2_possorted_genome_bam.bam"
  )

cell.list <- split(obj$seurat_clusters, obj$orig.ident)
str(cell.list)
List of 4
 $ sec1_ant: Factor w/ 17 levels "0","1","2","3",..: 13 13 2 2 2 2 2 13 2 2 ...
  ..- attr(*, "names")= chr [1:2695] "sec1_ant_TTGCCGGTGATCCCTC-1" "sec1_ant_CGGTATAGGTATTAGC-1" "sec1_ant_CGTTGAGTAATTGCGT-1" "sec1_ant_TAAATGCCGTCTCATG-1" ...
 $ sec1_pos: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 13 2 13 2 13 2 ...
  ..- attr(*, "names")= chr [1:3342] "sec1_pos_TGGGACCATTGGGAGT-1" "sec1_pos_CCGGTGCGAGTGATAG-1" "sec1_pos_TAGCCAGAGGGTCCGG-1" "sec1_pos_GAAGGGCATAACCATG-1" ...
 $ sec2_ant: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 2 2 13 2 13 2 ...
  ..- attr(*, "names")= chr [1:2823] "sec2_ant_TTGAGGCATTTAACTC-1" "sec2_ant_CACCGGAGATATCTCC-1" "sec2_ant_CACGTCTATGATGTGG-1" "sec2_ant_TGACCAGCTTCAAAGT-1" ...
 $ sec2_pos: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 2 2 2 2 2 2 ...
  ..- attr(*, "names")= chr [1:3288] "sec2_pos_TCCGCAGCCACCTAGC-1" "sec2_pos_GGTCCTTCATACGACT-1" "sec2_pos_GTGCCTAGCTATGCTT-1" "sec2_pos_GCGAGAAACGGGAGTT-1" ...
cell.list1 <- lapply(cell.list, function(x) {
  nm <- names(x)
  names(x) <- gsub(".*_(.*)","\\1",nm)
  x
})
str(cell.list1)
List of 4
 $ sec1_ant: Factor w/ 17 levels "0","1","2","3",..: 13 13 2 2 2 2 2 13 2 2 ...
  ..- attr(*, "names")= chr [1:2695] "TTGCCGGTGATCCCTC-1" "CGGTATAGGTATTAGC-1" "CGTTGAGTAATTGCGT-1" "TAAATGCCGTCTCATG-1" ...
 $ sec1_pos: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 13 2 13 2 13 2 ...
  ..- attr(*, "names")= chr [1:3342] "TGGGACCATTGGGAGT-1" "CCGGTGCGAGTGATAG-1" "TAGCCAGAGGGTCCGG-1" "GAAGGGCATAACCATG-1" ...
 $ sec2_ant: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 2 2 13 2 13 2 ...
  ..- attr(*, "names")= chr [1:2823] "TTGAGGCATTTAACTC-1" "CACCGGAGATATCTCC-1" "CACGTCTATGATGTGG-1" "TGACCAGCTTCAAAGT-1" ...
 $ sec2_pos: Factor w/ 17 levels "0","1","2","3",..: 13 13 13 2 2 2 2 2 2 2 ...
  ..- attr(*, "names")= chr [1:3288] "TCCGCAGCCACCTAGC-1" "GGTCCTTCATACGACT-1" "GTGCCTAGCTATGCTT-1" "GCGAGAAACGGGAGTT-1" ...
TrackPlot(bamfile = bamfiles, cell.group = cell.list1, gtf = gtf, gene = "Sept8", highlights = c(53548291,53549565), junc=TRUE)

Back to top