Package 'seqpac'

Description

add_reanno adds biotypes from reannotation to PAC.

Usage

add_reanno(
  reanno,
  mismatches = 0,
  type = "genome",
  bio_search,
  bio_perfect = FALSE,
  genome_max = 10,
  merge_pac = NULL
)

Arguments

Details

Given a reanno object generated by make_reanno this function will either extract genome coordinates or classify biotypes using a list of search terms (see details). Information will be compiled into a data frame with the same order of sequences as in the original PAC master file.

Please rename the column names of your previously added columns if running map_reanno more than two times using the same type input ("genome" or "biotype").

Value

Data frame with mismatch and coordinate or biotype information (search term hits) that can directly be added to an Anno data.frame of a PAC object, given that the PAC object was used in the reannotation workflow. Annotations can be further simplified using the simplify_reanno function.

See Also

http://bowtie-bio.sourceforge.net/index.shtml for information about Bowtie and https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: PAC_annotate(), as.reanno(), import_reanno(), make_conv(), make_reanno(), map_reanno(), simplify_reanno()

Examples

######################################################### 
##### Example for reference mapping and classification 
closeAllConnections()
######################################################### 
##### Create an reanno object

### First, if you haven't already generated Bowtie indexes for the included
# fasta references you need to do so. If you are having problem see the small
# RNA guide (vignette) for more info.
                                                             
 ## tRNA:
 trna_file <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 trna_dir <- dirname(trna_file)
 
 if(!sum(stringr::str_count(list.files(trna_dir), ".ebwt")) ==6){     
     Rbowtie::bowtie_build(trna_file, 
                           outdir=trna_dir, 
                           prefix="tRNA", force=TRUE)
                           }
 ## rRNA:
 rrna_file <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 rrna_dir <- dirname(rrna_file)
 
 if(!sum(stringr::str_count(list.files(rrna_dir), ".ebwt")) ==6){
     Rbowtie::bowtie_build(rrna_file, 
                           outdir=rrna_dir, 
                           prefix="rRNA", force=TRUE)
                           }
                           
                           
##  Then load a PAC-object and remove previous mapping from anno:
 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 ref_paths <- list(trna= trna_file, rrna= rrna_file)

                                    
##  You may add an output path of your choice, but here we use a temp folder:
 output <- paste0(tempdir(),"/seqpac/test")


##  Then map the PAC-object against the fasta references. 
# Warning: if you use your own data, you may want to use override=FALSE, to avoid
# deleting previous mapping by mistake.

map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=0,  import="biotype", 
               threads=2, keep_temp=FALSE, override=TRUE)
 
##  Then import and generate a reanno-object of the temporary bowtie-files
reanno_biotype <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
 
                                    
## Now make some search terms against reference names to create shorter names
# Theses can be used to create factors in downstream analysis
# One search hit (regular expressions) gives one new short name 

bio_search <- list(
              rrna=c("5S", "5.8S", "12S", "16S", "18S", "28S", "pre_S45"),
              trna =c("_tRNA", "mt:tRNA"))
 
 
## You can merge directly with your PAC-object by adding your original 
# PAC-object, that you used with map_reanno, to merge_pac option.

pac <- add_reanno(reanno_biotype, bio_search=bio_search, 
                       type="biotype", bio_perfect=FALSE, 
                       mismatches = 0, merge_pac=pac)

table(anno(pac)$mis0_bio) 

############################################################################ 
## Similarly, you can use Bowtie indexed genome fasta references
## But don't forget to use import="genome" for coordinate import
#   
# ref_paths <- list(genome="<path_to_bowtie_indexed_fasta>")
# output_genome <- "<your_path_to_output_folder>"
#
## We can actually map against several genomes simultaneously:
# (to illustrate the principle we run the mycoplasma genome twice)
# ref_paths <- list(genome1=mycoplasma_file, genome2=mycoplasma_file)
# ### Run map_reanno 
# map_reanno(PAC=pac, input=ref_paths, output=output, 
#            type="internal", mismatches=0, import="genome", 
#            threads=2, keep_temp=TRUE, override=TRUE)
# reanno_genome <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
############################################################################
#### The trick to succeed with bio_perfect=TRUE 

## Run add_reanno with bio_perfect="FALSE" (look where "Other=XX" occurs)
# Use the bio
#' ## Add a search terms that catches the other rrna 

anno <- add_reanno(reanno_biotype, bio_search=bio_search, 
                   type="biotype", bio_perfect=FALSE, mismatches = 0)

## Find sequences that has been classified as other
other_seqs  <- anno[grepl("other", anno$mis0),]$seq
tab <- full(reanno_biotype)$mis0$trna
tab[tab$seq %in% other_seqs,]         #No other hit in trna

tab <- full(reanno_biotype)$mis0$rrna
tab[tab$seq %in% other_seqs,]         #A few in rrna


## Add a search terms that catches the other rrna 
bio_search <- list(
              rrna=c("5S", "5.8S", "12S", "16S", 
                     "18S", "28S", "pre_S45", "Other"),
              trna =c("_tRNA", "mt:tRNA"))
                
anno <- add_reanno(reanno_biotype, bio_search=bio_search, 
                   type="biotype", bio_perfect=FALSE, mismatches = 0)

table(anno$mis0)                 

## Repeat search until no "Other" appear when running add_reanno, 
## then run  bio_perfect=TRUE: 

anno <- add_reanno(reanno_biotype, bio_search=bio_search,  
                   type="biotype", bio_perfect=TRUE, mismatches = 0)

Description

anno(PAC)

Access the annotation table in S4 PAC

Usage

anno(object)

## S4 method for signature 'PAC'
anno(object)

## S4 replacement method for signature 'PAC,data.frame'
anno(object) <- value

Arguments

Value

The Anno data.frame table from a PAC-object.

Updated S4 PAC-object.

See Also

Other PAC methods: counts(), names,PAC-method, norm(), pheno(), show.PAC(), summary()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

as.PAC Converts an S3 PAC object (list) to an S4 PAC object

Usage

as.PAC(from)

Arguments

Details

Seqpac comes with two versions of the PAC object, either S4 or S3. The S3 PAC is simply a list where each object can be received using the $-sign. The S4 version is a newer type of R object with predefined slots, that is received using the @-sign. This function converts an S3 PAC object (list) to an S4 PAC object. You can also use the S4 coercion method "as" to turn an S4 PAC into an S3. See examples below.

Value

An S4 PAC object.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))
                  
# check type
class(pac)
isS4(pac)
 
 
# Turns S4 PAC object into a S3
pac_s3 <- as(pac, "list")

# Turns S3 PAC object into a S4
pac_s4 <- as.PAC(pac_s3)

Description

as.reanno Converts an S3 reanno object (list) to an S4 reanno object

Usage

as.reanno(from)

Arguments

Details

Seqpac comes with two versions of the reanno object, either S4 or S3. The S3 is simply a list where each object can be received using the $-sign. The S4 version is a newer type of R object with predefined slots, that is received using the @-sign. This function converts an S3 reanno object (list) to an S4 object. You can also use the S4 coercion method "as" to turn an S4 into an S3. See examples below.

Value

An S4 reanno object.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: PAC_annotate(), add_reanno(), import_reanno(), make_conv(), make_reanno(), map_reanno(), simplify_reanno()

Examples

##### Create an reanno object

##  First load a PAC- object

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 
##  Then specify paths to fasta references
# If you are having problem see the vignette small RNA guide for more info.
 
 trna_file <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 trna_dir <- dirname(trna_file)
 
 if(!sum(stringr::str_count(list.files(trna_dir), ".ebwt")) ==6){     
     Rbowtie::bowtie_build(trna_file, 
                           outdir=trna_dir, 
                           prefix="tRNA", force=TRUE)
                           }
 
 ref_paths <- list(trna= trna_file)
                                    
##  Add output path of your choice.
# Here we use the R temporary folder depending on platform
 output <- paste0(tempdir(),"/seqpac/test")

## Make sure it is empty (otherwise you will be prompted for a question)
out_fls  <- list.files(output, recursive=TRUE)
closeAllConnections()
suppressWarnings(file.remove(paste(output, out_fls, sep="/")))

##  Then map your PAC-object against the fasta references
 map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=0,  import="biotype", 
               threads=2, keep_temp=FALSE, override=TRUE)

##  Then import and generate a reanno-object of the temporary bowtie-files
reanno <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
                                    
## Now make some search terms against reference names to create shorter names
# Theses can be used to create factors in downstream analysis
# One search hit (regular expressions) gives one new short name 
bio_search <- list(trna =c("_tRNA", "mt:tRNA"))
 
 
## You can merge directly with your PAC-object by adding your original 
# PAC-object, that you used with map_reanno, to merge_pac option.
 pac <- add_reanno(reanno, bio_search=bio_search, 
                       type="biotype", bio_perfect=FALSE, 
                       mismatches = 0, merge_pac=pac)
                       
                       
## Turn your S3 list to an S4 reanno-object
isS4(reanno)
names(reanno)
table(overview(reanno)$trna)
 
reanno_s3 <- as(reanno, "list")
isS4(reanno_s3)   

# Turns S3 reanno object into a S4 
reanno_s4 <- as.reanno(reanno_s3)
isS4(reanno_s4)

Description

sample names

Usage

## S4 method for signature 'PAC'
colnames(x)

Arguments

Value

The sample names held by the Counts table (column names) and the Pheno table (row names) of a PAC-object.

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

Access the raw counts table in S4 PAC

Usage

counts(object)

## S4 method for signature 'PAC'
counts(object)

## S4 replacement method for signature 'PAC,data.frame'
counts(object) <- value

Arguments

Value

The Counts data frame table from a PAC-object.

Updated S4 PAC-object.

See Also

Other PAC methods: anno(), names,PAC-method, norm(), pheno(), show.PAC(), summary()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

filtsep_bin Converts PAC_filtsep data.frame output into a binary (hit-no-hit) data.frame

Usage

filtsep_bin(filtsep_out)

Arguments

Details

Given a PAC_filtsep output data.frame, where each column contains the sequences that passed the filter for a specific group specified in pheno_target, filtsep_bin converts this into a data.frame where sequences are reported as hit (=1) or no hit (=0). Such binary coded group occurrence can for example be used by UpSetR::upset to generate visualization of overlaps using UpSet plots.

Value

data.frame where each uniques sequence (row names) in filtsep_out are reported as hit (=1) or no hit (=0) across samples or sample groups (column names).

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), map_rangetype(), tRNA_class()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

## Keep sequences with 5 counts (threshold) in 100% (coverage) of 
## samples in a group:
 # Use PAC_filtsep to find sequences 
 filtsep <- PAC_filtsep(pac, norm="counts", threshold=5, 
                        coverage=100, pheno_target= list("stage"))
                        
 # Filter by unique sequences passing filtsep  
 filtsep <- unique(do.call("c", as.list(filtsep)))
 pac_filt <- PAC_filter(pac, subset_only = TRUE, anno_target= filtsep)
 
 # Find overlap
 olap <- reshape2::melt(filtsep, 
                        measure.vars = c("Stage1", "Stage3", "Stage5"), 
                        na.rm=TRUE)
                        
## Upset plot using the UpSetR package
 # (when output="binary" PAC_filtsep uses filtsep_bin for binary conversion
 # Use PAC_filtsep with binary output
 filtsep_bin <- PAC_filtsep(pac, norm="counts", threshold=5, 
                            coverage=100, pheno_target= list("stage"), 
                            output="binary")
 
# Plot Venn diagram or UpSetR
#
# plot(venneuler::venneuler(data.frame(olap[,2], olap[,1]))) 
#
# UpSetR::upset(filtsep_bin, sets = colnames(filtsep_bin), 
#              mb.ratio = c(0.55, 0.45), order.by = "freq", keep.order=TRUE)

Description

Access the multilevel list of imported data.frames in a S4 reanno object

Usage

full(object)

## S4 method for signature 'reanno'
full(object)

## S4 replacement method for signature 'reanno,list'
full(object) <- value

Arguments

Value

All the results that were imported into the reanno-object returned as a list of tibble data.frames (class tbl_df, tbl, data.frame).

Updated S4 reanno-object.

See Also

Other reanno methods: overview()

Examples

######################################################### 
##### Create a reanno object

##  First load a PAC- object

 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 
##  Then specify paths to fasta references
# If you are having problem see the vignette small RNA guide for more info.
 
 trna_path <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)  
 rrna_path <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 
 ref_paths <- list(trna= trna_path, rrna= rrna_path)
 
##  Add output path of your choice.
# Here we use the R temporary folder depending on platform  
if(grepl("windows", .Platform$OS.type)){
 output <- paste0(tempdir(), "\\seqpac\\test")
}else{
 output <- paste0(tempdir(), "/seqpac/test")}

## Make sure it is empty (otherwise you will be prompted for a question)
out_fls  <- list.files(output, recursive=TRUE)
closeAllConnections()
suppressWarnings(file.remove(paste(output, out_fls, sep="/")))

##  Then map your PAC-object against the fasta references 
 map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=2,  import="biotype", 
               threads=2, keep_temp=FALSE)

##  Then generate a reanno-object of the temporary bowtie-files 
reanno_object <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)

## Accessing content and S4/S3 conversion:
names(reanno_object)
overview(reanno_object)
full(reanno_object)
rownames(reanno_object)
length(reanno_object)
nrow(reanno_object)
reanno_s3 <- as(reanno_object, "list")
reanno_s4 <- as.reanno(reanno_s3)

Description

This function imports bowtie output and summarizes the content.

Usage

import_reanno(
  bowtie_path,
  threads = 1,
  coord = FALSE,
  report = "minimum",
  reduce = NULL
)

Arguments

Details

Given the path to alignment outputs from bowtie, import_reanno will attempt to read these files into R and generate a list of unsorted data.frames where each row summarizes all annotations for a given sequence. It is called by map_reanno to generate summarized data from bowtie output saved as .Rdata files. These files can then be re-imported into R where information is organized using make_reanno and finally added to the annotation table (PAC$Anno) of the original PAC-list object using add_reanno.

Value

List of data frames with additional information from reannotation files generated with bowtie. If report="minimum", the function will report each hit only as the number of mismatches for a given reference file. If report="full" the full name reported in the fasta file used as reference in the bowtie reannotation will be reported. If a reference name is specified in reduce, this reference is excerpted from report="full" and is instead reported as report="minimum".

Caution: Large references with lots of redundancy (such as pirBase in some species) will cause massive character strings if report="full" with no restrictions. Specifying such references in reduce=<reference_names> will circumvent this problem.

Important

Re-annotation must have been done using Bowtie default output settings, where output files specifically have been named '.txt'. The function will not work with other formats such as SAM/BAM formats.

The basenames of the bowtie output will also be used to annotate. Example: bowtie -v 3 -a -f '<piRBase.fasta>' '<master_anno.fasta>' piRNA.txt Will annotate sequences as piRNA since the txt-file output was named 'piRNA'

See Also

http://bowtie-bio.sourceforge.net/index.shtml for information about Bowtie and for Rbowtie: https://www.bioconductor.org/packages/release/bioc/html/Rbowtie.html. https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: PAC_annotate(), add_reanno(), as.reanno(), make_conv(), make_reanno(), map_reanno(), simplify_reanno()

Examples

######################################################### 
##### Test import_reanno

### First, if you haven't already generated Bowtie indexes for the included
# fasta reference you need to do so. If you are having problem see the small
# RNA guide (vignette) for more info.

## tRNA reference:
trna_file <- system.file("extdata/trna", "tRNA.fa", 
                         package = "seqpac", mustWork = TRUE)
trna_dir<- gsub("tRNA.fa", "", trna_file)

if(!sum(stringr::str_count(list.files(trna_dir), ".ebwt")) ==6){
  Rbowtie::bowtie_build(trna_file, 
                        outdir=trna_dir, 
                        prefix="tRNA", force=TRUE)
}
##  Then load a PAC-object and remove previous mapping from anno:
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                 package = "seqpac", mustWork = TRUE))

ref_paths <- list(trna= trna_file)

##  You may add an output path of your choice, but here we use a temp folder:
output <- paste0(tempdir(),"/seqpac/test")

##  Then map the PAC-object against the fasta references. Warning: if you use
# your own data, you may want to use override=FALSE, to avoid deleting previous
# mapping by mistake. keep_temp=TRUE can be used to run import_reanno
# independently.

map_reanno(pac, input=ref_paths, output=output,
           type="internal", mismatches=2,  import="biotype", 
           threads=2, keep_temp=TRUE, override=TRUE)

reanno1 <- import_reanno(output, report="Biotype",  threads=1)

Description

number of objects in PAC

length of reanno object

Usage

## S4 method for signature 'PAC'
length(x)

## S4 method for signature 'reanno'
length(x)

Arguments

Value

The number of items in a PAC-object.

Number of items in the reanno-object.

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra functionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

Generates a table that can be used for converting seqnames across referens genomes.

Usage

make_conv(
  reference_list = NULL,
  ref_path_A = NULL,
  ref_path_B = NULL,
  ref_path_C = NULL,
  output = "tibble",
  skip_after = NULL
)

Arguments

Details

Given the file paths to 2-3 fasta reference genomes, this function will quickly match the the sequences between these reference genomes and return seqnames of each matching sequence. Thus, given the paths to the fasta for the Ensembl, UCSC and NCBI assemblies of a given genome version (e.g. hg38 or dm6), will result in a name conversion table between these databases.

Only perfect matches or no matches will be reported. Thus, in case an entire sequence is missing (e.g. a sex chromosomes), a table will be returned with a warning, but if a sequence is partly missing (e.g. only half a chromosome), then an error is returned.

Sequence matching is done using md5 hashes, which dramatically increases the speed for perfect matches. Matching will always be done from reference A against the other references. Thus, if reference A (=1st reference in reference_list) contains less sequences than the other references, only the reference A sequences will be reported in the output, having the same order as in reference A.

Value

A name conversion table either as a tibble or a data frame (see output)

See Also

https://github.com/OestLab/seqpac for updates.

Other PAC reannotation: PAC_annotate(), add_reanno(), as.reanno(), import_reanno(), make_reanno(), map_reanno(), simplify_reanno()

Examples

## Only for testing:
fasta_path <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE) 
ref1 <- Biostrings::readDNAStringSet(fasta_path) 
ref1 <- ref1[1:295]
sqnames <- do.call("rbind",(strsplit(names(ref1), "\\d chr")))[,2]
names(ref1) <- do.call("rbind",(strsplit(sqnames, " \\(")))[,1]
logi_dup <- duplicated(do.call("rbind", strsplit(names(ref1),"\\:"))[,1])
ref1 <- ref1[!logi_dup]
ref2 <- ref1
names(ref2) <- paste0("chr", names(ref2))
# Save new reference in temporary folder 
if(grepl("windows", .Platform$OS.type)){
  tmpdr <- paste0(tempdir(), "\\seqpac")
}else{
  tmpdr <- paste0(tempdir(), "/seqpac")}
dir.create(tmpdr, showWarnings=FALSE) 
ref_path1 <- paste0(tmpdr, "/ref1.fa")
ref_path2 <- paste0(tmpdr, "/ref2.fa") 
Biostrings::writeXStringSet(ref1, filepath=ref_path1, format="fasta")
Biostrings::writeXStringSet(ref2, filepath=ref_path2, format="fasta")
ref_list <- list(ensembl=ref_path1, ucsc=ref_path2)
conv_table <- make_conv(reference_list=ref_list)
conv_table
              

## The principles:
#
#ref_path_A <- "/some/path/to/ensembl.fa"
#ref_path_B <- "/some/path/to/ucsc.fa"
#ref_path_C <- "/some/path/to/refseq.fa"
#
#ref_list <- list(ensembl=ref_path_A, UCSC=ref_path_B, NCBI=ref_path_C)
#
## Best (user defined names):
#conv_table <- make_conv(reference_list=ref_list) 
#
## But also (no names)
#conv_table <- make_conv(ref_path_A, ref_path_B, ref_path_C)
#conv_table <- make_conv(ref_path_A, ref_path_C)
#
#
## Make short names (skip everything after white space)
#conv_table <- make_conv(reference_list=reference_list, skip_after=" ")

Description

make_counts uses parallel processing to generate a count table.

Usage

make_counts(
  input,
  trimming = NULL,
  threads = 1,
  save_temp = FALSE,
  plot = TRUE,
  parse = "default_illumina",
  evidence = c(experiment = 2, sample = 1),
  optimize = list(on_disk = FALSE, chunk_size = NULL, rm_wild = FALSE, seq_min = NULL)
)

Arguments

Details

Given a path to fastq files this function performs low-level evidence filtering, generating a counts table of sequences passing the filter and plots summary statistics.

Value

A list containing three objects:

1. counts (data frame) = count table.

2. progress_report = progress report from trimming and evidence filter.

3. evidence_plots = bar graphs showing the impact of evidence filter.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), PAC_create(), make_PAC(), make_cutadapt(), make_pheno(), make_trim(), merge_lanes()

Examples

############################################################ 
### Seqpac fastq trimming with the make_counts function 
### using default settings for NEBNext small RNA adaptor 

# Seqpac includes strongly down-sampled smallRNA fastq.
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
input <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)

# Now we can run make_counts
# Notice that make_counts will generate another temp folder, that will 
# be emptied on finalization. By setting save_temp=TRUE you may save the 
# content. You may also use your own trimming settings by using the parse 
# option. See ?make_trim for more details. 
 
counts  <- make_counts(input, threads=1, parse="default_neb",
                       trimming="seqpac", plot=TRUE,
                       evidence=c(experiment=2, sample=1))

head(counts$counts)
head(counts$progress_report)

# Notice that that there are fewer unique sequences than number of reads 
# passing the filter. In normal, large, fastq we keep 95-99% of all reads 
# while filtering out 30-40% of the sequences that only appear in 1 fastq 
# Thus, the evidence filter may gain you performance in later analysis by 
# avoiding nonsense sequences. 
     
 
 #############################################################
 ### Lets change the evidence filter
 ###
 
 # 2 evidence over two indepenent samples, saving single sample 
 # sequences reaching 3 counts:
 
 test <- make_counts(input=input, trimming="seqpac", 
                     parse="default_neb",  
                     evidence=c(experiment=2, sample=3))
 extras <- apply(test$counts, 1, function(x){sum(!x==0)})
 test$counts[extras==1,]  # 6 single sample sequences reached 3 counts
 
 # 2 evidence over two independent samples, saving single sample 
 # sequences reaching 2 counts
  
 test <- make_counts(input=input,  trimming="seqpac", 
                     parse="default_neb",  
                     evidence=c(experiment=2, sample=2))
 extras <- apply(test$counts, 1, function(x){sum(!x==0)})
 test$counts[extras==1,] # 120 single sample sequences reached 2 counts

Description

make_cuadapt cutadapt/fastq_quality_filter

Usage

make_cutadapt(input, output, parse = NULL, threads = 1)

Arguments

Details

Given a path to sequence files in fastq format this function will trim adaptor and remove sequences with low quality.

Value

Externally the function will generate trimmed and/or quality filtered fastq files in the output folder. Internally, a list of logs that can be used to generate a progress report is returned.

See Also

https://cutadapt.readthedocs.io/en/stable/ for download and documentation on cutadapt. http://hannonlab.cshl.edu/fastx_toolkit/commandline.html for download and documentation on fastq_quality_filter. https://github.com/OestLab/seqpac for updates on seqpac.

Other PAC generation: PAC_check(), PAC_create(), make_PAC(), make_counts(), make_pheno(), make_trim(), merge_lanes()

Examples

############################################################      
### Principle of trimming using the make_cutadapt function
### (Important: Need external installations of cutadapt 
###  and fastq_quality_filter to work)
#  
#   input = "/some/path/to/input/folder"
#   output =  "/some/path/to/output/folder"
# 
## Parse for make_cutadapt is a list of 2 character string expressions.
## The first is parsed to cutadapt and the other to fastq_quality_filter 
## For parallel processes '-j 1' is recommended since seqpac will   
## parallelize across samples and not within.
## Run system2("cutadapt -h", stdout=TRUE) and 
## system("fastq_quality_filter -h", stdout=TRUE) 
## for more options.
#   
## String to parse to cutadapt:
# cut_prs <- paste0("-j 1 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCACAT",
#                    " --discard-untrimmed --nextseq-trim=20",
#                    " -O 10 -m 7 -M 70")
# 
## Add string to parse to fastq_quality_filter:
#  parse = list(
#            cutadapt=cut_prs,
#            fastq_quality_filter="-q 20 -p 80")
#               
#  logs  <-  make_cutadapt(input, output, threads=8, parse=parse)

#' # Clean up temp
closeAllConnections()
fls_temp  <- list.files(tempdir(), recursive=TRUE, full.names = TRUE)
file.remove(fls_temp, showWarnings=FALSE)

Description

make_PAC Compiles, order and checks the assumptions of a PAC file.

Usage

make_PAC(counts, pheno = NULL, anno = NULL, output = "S4")

Arguments

Details

Given a pheno, an anno and a counts data.frames the functions will safely generate a simple PAC list-object.

Value

Ordered PAC list object, checked for compatibility with downstream analysis in seqpac.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), PAC_create(), make_counts(), make_cutadapt(), make_pheno(), make_trim(), merge_lanes()

Examples

### First make counts 

# Seqpac includes strongly down-sampled smallRNA fastq.
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
input <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)

# Notice that make_counts will generate another temp folder, that will 
# be emptied on finalization. By setting save_temp=TRUE you may save the 
# content.  
 
counts  <- make_counts(input, threads=2, parse="default_neb",
                       trimming="seqpac", plot=TRUE,
                       evidence=c(experiment=2, sample=1))

colnames(counts$counts)

### Then generate a phenotype table with make_pheno

#  Note:  'Sample_ID' column needs to be similar IDs as 
#          colnames in the counts table. You may also 
#          specify a path to a txt file.

Sample_ID <- colnames(counts$counts)

pheno <- data.frame(Sample_ID=Sample_ID,
                       Treatment=c(rep("heat", times=1), 
                                   rep("control", times=2)),
                       Batch=rep(c("1", "2", "3"), times=1)) 

pheno <- make_pheno(pheno=pheno, progress_report=counts$progress_report, 
                     counts=counts$counts)


pheno 

# Note that progress report from make_counts is added if you specify it  
     
### Lastly combine into PAC

pac <- make_PAC(pheno=pheno, counts=counts$counts)


pac
names(pac)

# Note: a simple annotation table is added automatically.
head(anno(pac))

# Clean up temp
closeAllConnections()
fls_temp  <- list.files(tempdir(), recursive=TRUE, full.names = TRUE)
file.remove(fls_temp, showWarnings=FALSE)

Description

make_pheno generates a Pheno data.frame object containing additional metadata for unique samples.

Usage

make_pheno(pheno, counts = NULL, progress_report = NULL)

Arguments

Details

Given the path to a directory with a single csv file, the function will read the file. Alternatively, a data.frame with sample information can be provided. If provided with a count table, (see make_counts), it will attempt to organize the row names according to the column names in the counts table.

Value

A Pheno data.frame compatible with make_PAC

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), PAC_create(), make_PAC(), make_counts(), make_cutadapt(), make_trim(), merge_lanes()

Examples

### First make counts 

# Seqpac includes strongly down-sampled smallRNA fastq.
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
input <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)

# Notice that make_counts will generate another temp folder, that will 
# be emptied on finalization. By setting save_temp=TRUE you may save the 
# content.  
 
counts  <- make_counts(input, threads=2, parse="default_neb",
                       trimming="seqpac", plot=TRUE,
                       evidence=c(experiment=2, sample=1))

colnames(counts$counts)

### Then generate a phenotype table with make_pheno

#  Note:  'Sample_ID' column needs to be similar IDs as 
#          colnames in the counts table. You may also 
#          specify a path to a txt file.

Sample_ID <- colnames(counts$counts)

pheno <- data.frame(Sample_ID=Sample_ID,
                       Treatment=c(rep("heat", times=1), 
                                   rep("control", times=2)),
                       Batch=rep(c("1", "2", "3"), times=1)) 

pheno <- make_pheno(pheno=pheno, progress_report=counts$progress_report, 
                     counts=counts$counts)


pheno 

# Note that progress report from make_counts is added if you specify it  
     
### Lastly combine into PAC

pac <- make_PAC(pheno=pheno, counts=counts$counts)


pac
names(pac)

# Note: a simple annotation table is added automatically.
head(anno(pac))

# Clean up temp
closeAllConnections()
fls_temp  <- list.files(tempdir(), recursive=TRUE, full.names = TRUE)
file.remove(fls_temp, showWarnings=FALSE)

Description

make_reanno makes a reannotation (reanno) object.

Usage

make_reanno(reanno_path, PAC, mis_fasta_check = FALSE, output = "S4")

Arguments

Details

Given the path to the R reannotation files (reanno_mis0/1/2/3.Rdata) generated by map_reanno, this function will summarize the reannotation files into one output that matches the order of sequences in a PAC object.

Value

Contains two slots/objects: 1. Overview: A table with a summary of the mapping. 2. Full_anno: Lists of tables holding the full imports per mismatch cycle (mis0, mis1 etc) and reference (mi0-ref1, mis0-ref2, mis1-ref1, mis1-ref2 etc). All table rows are ordered according to the order of the PAC originally used for the mapping. If mis_fasta_check is specified, the function will look for a anno_misX.fa (X = the file with the highest number) previously generated by the reannotation workflow. This file is used to double check so that no sequences are missing before and after reannotation.

See Also

http://bowtie-bio.sourceforge.net/index.shtml for information about Bowtie https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: PAC_annotate(), add_reanno(), as.reanno(), import_reanno(), make_conv(), map_reanno(), simplify_reanno()

Examples

######################################################### 
##### Simple example for reference mapping 
##### Please see manual for ?simply_reanno an ?add_reanno for more advanced 
##### mapping
closeAllConnections()
######################################################### 
##### Create an reanno object

### First, if you haven't already generated Bowtie indexes for the included
# fasta references you need to do so. If you are having problem see the small
# RNA guide (vignette) for more info.
                                                             
 ## tRNA:
 trna_file <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 trna_dir <- dirname(trna_file)
 
 if(!sum(stringr::str_count(list.files(trna_dir), ".ebwt")) ==6){     
     Rbowtie::bowtie_build(trna_file, 
                           outdir=trna_dir, 
                           prefix="tRNA", force=TRUE)
                           }
 ## rRNA:
 rrna_file <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 rrna_dir <- dirname(rrna_file)
 
 if(!sum(stringr::str_count(list.files(rrna_dir), ".ebwt")) ==6){
     Rbowtie::bowtie_build(rrna_file, 
                           outdir=rrna_dir, 
                           prefix="rRNA", force=TRUE)
                           }
                           
                           
##  Then load a PAC-object and remove previous mapping from anno:
 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 ref_paths <- list(trna= trna_file, rrna= rrna_file)

                                    
##  You may add an output path of your choice, but here we use a temp folder:
 output <- paste0(tempdir(),"/seqpac/test")


##  Then map the PAC-object against the fasta references. 
# Warning: if you use your own data, you may want to use override=FALSE, to avoid
# deleting previous mapping by mistake.

map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=0,  import="biotype", 
               threads=2, keep_temp=FALSE, override=TRUE)
 
##  Then import and generate a reanno-object of the temporary bowtie-files
reanno_biotype <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
 
                                    
## Now make some search terms against reference names to create shorter names
# Theses can be used to create factors in downstream analysis
# One search hit (regular expressions) gives one new short name 

bio_search <- list(
              rrna=c("5S", "5.8S", "12S", "16S", "18S", "28S", "pre_S45"),
              trna =c("_tRNA", "mt:tRNA"))
 
 
## You can merge directly with your PAC-object by adding your original 
# PAC-object, that you used with map_reanno, to merge_pac option.

pac <- add_reanno(reanno_biotype, bio_search=bio_search, 
                       type="biotype", bio_perfect=FALSE, 
                       mismatches = 0, merge_pac=pac)

table(anno(pac)$mis0_bio) 

############################################################################ 
## Similarly, you can use Bowtie indexed genome fasta references
## But don't forget to use import="genome" for coordinate import
#    
# ref_paths <- list(genome="<path_to_bowtie_indexed_fasta>")
# output_genome <- "<your_path_to_output_folder>"

## We can actually map against several genomes simultaneously:
# (to illustrate the principle we run the mycoplasma genome twice)

 mycoplasma_file <- system.file("extdata/mycoplasma_genome",
                                "mycoplasma.fa",
                                package = "seqpac", mustWork = TRUE)
 mycoplasma_dir<- gsub("mycoplasma.fa", "", mycoplasma_file)
 
## Don't forget to create bowtie indexes
 if(!sum(stringr::str_count(list.files(mycoplasma_dir), ".ebwt")) ==6){
     Rbowtie::bowtie_build(mycoplasma_file, 
                           outdir=mycoplasma_dir, 
                           prefix="mycoplasma", force=TRUE)
                           }
                           
 ref_paths <- list(genome1=mycoplasma_file, genome2=mycoplasma_file)
 map_reanno(PAC=pac, input=ref_paths, output=output, 
            type="internal", mismatches=0, import="genome", 
            threads=2, keep_temp=TRUE, override=TRUE)
 reanno_genome <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
 
 table(overview(reanno_genome)$Any) # Low Mycoplasma contamination
 
############################################################################

Description

make_trim uses parallel processing to generate trimmed reads.

Usage

make_trim(
  input,
  output,
  indels = TRUE,
  concat = 12,
  check_mem = FALSE,
  threads = 1,
  chunk_size = NULL,
  polyG = c(type = NULL, min = NULL, mismatch = NULL),
  adapt_3_set = c(type = "hard_rm", min = 10, mismatch = 0.1),
  adapt_3 = "AGATCGGAAGAGCACACGTCTGAACTCCAGTCACTA",
  adapt_5_set = c(type = NULL, min = NULL, mismatch = NULL),
  adapt_5 = NULL,
  nucleotide_range = c(min = NULL, max = NULL),
  quality = c(threshold = 20, percent = 0.8)
)

Arguments

Details

Given a file path to sequence files in fastq format (extension .fastq.gz or .fastq) this function will first trim the sequences from adaptor sequence, and then apply sequence size and phred score quality filters. On systems with low memory or very large fastq files, users may choose to handle each file in chunks on the expense of time performance.

Value

Exports a compressed trimmed fastq file to the output folder (file name: basename+trim.fastq.gz). If any type="save" has been set, an additional fastq file (file name includes 'REMOVED') will be exported. In addition, overview statisics (progress report) will be returned as a data.frame.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), PAC_create(), make_PAC(), make_counts(), make_cutadapt(), make_pheno(), merge_lanes()

Examples

############################################################ 
### Seqpac fastq trimming with the make_trim function 
### (for more streamline options see the make_counts function) 
 ############################################################ 
### Seqpac fastq trimming with the make_counts function 
### using default settings for NEBNext small RNA adaptor 

# Seqpac includes strongly down-sampled smallRNA fastq.
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
input <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)

# Run make_trim using NEB-next adaptor
# (Here we store the trimmed files in the input folder)
# (but you may choose where to store them using the output option)

input <- unique(dirname(input))  # mak_trim searches for fastq in input dir
output <- tempdir()              # save in tempdir


list.files(input, pattern = "fastq") #before

prog_report  <-  make_trim(
       input=input, output=output, threads=1, 
       adapt_3_set=c(type="hard_rm", min=10, mismatch=0.1), 
       adapt_3="AGATCGGAAGAGCACACGTCTGAACTCCAGTCACTA", 
       polyG=c(type="hard_trim", min=20, mismatch=0.1),
       nucleotide_range=c(min=14, max=70),
       quality=c(threshold=20, percent=0.8))
       
list.files(path = output, pattern = "fastq") #after

# How did it go? Check progress report:  
prog_report     
       

## The principle:
# input <- "/some/path/to/untrimmed_fastq_folder"
# output <- "/some/path/to/output_folder_for_trimmed_fastq"
# 
# prog_report  <-  make_trim(
#      input=input, output=output, 
#      threads=<how_many_parallel>, 
#      check_mem=<should_memory_be_checked?>, 
#      adapt_3_set=<options_for_main_trimming, 
#      adapt_3=<sequence_to be trimmed, 
#      polyG=<Illumina_Nextseq_type_poly_G_trimming>,
#      nucleotide_range=<what_sequence_range_to_save>,
#      quality=<quality_filtering_options>))

Description

map_rangetype Classifies sequences based on interval mapping against a reference.

Usage

map_rangetype(
  map,
  type = "percent",
  ss = NULL,
  min_loop_width = 4,
  intervals = list(start = 1:5, mid = 45:55, end = 95:100),
  N_include = FALSE
)

Arguments

Details

Given a PAC_map object (PAC_mapper) and an interval list this function will attempt to classify mapped sequences based on where these sequences starts and ends in reference. This function can for example be used for 5' and 3' tRNA classification.

Value

Map list object containing reference sequence (Ref_seq) as Biostrings::DNAStringSet and the new classifications embedded with the alignments (Alignments) in a dataframe.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), tRNA_class()

Examples

###########################################################
### test the map_rangetype function
# More complicated examples can be found in the vignette.
##----------------------------------------

# First create an annotation blank PAC with group means
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                 package = "seqpac", mustWork = TRUE))
anno(pac) <- anno(pac)[,1, drop=FALSE]
pac_trna <- PAC_summary(pac, norm = "cpm", type = "means", 
                        pheno_target=list("stage"), merge_pac = TRUE)

# Then re-annotate only tRNA using the PAC_mapper function
ref <- system.file("extdata/trna", "tRNA.fa", 
                         package = "seqpac", mustWork = TRUE)
map_object <- PAC_mapper(pac_trna, input=ref, N_up = "NNN", N_down = "NNN", 
                         mismatches=0, threads=2, report_string=TRUE, 
                         override=TRUE)

## Coverage plot of tRNA using PAC_covplot

# Single tRNA targeting a summary table 
PAC_covplot(pac_trna, map=map_object, 
                      summary_target= list("cpmMeans_stage"),
                      map_target="tRNA-Ala-AGC-1-1")
            
## Classify range types with map_rangetype (see vignette for examples
# on how to use ss-files for detailed tRNA loop structure).

# Classify fragments using percent intervals
map_object <- map_rangetype(map_object, 
                intervals = list(start = 1:5, mid = 45:55, end = 95:100))
       
names(map_object)
map_object[[1]]

Description

Aligns sequences in a PAC list object against specified references generating summarized report files in the destination folder.

Usage

map_reanno(
  PAC,
  type = "internal",
  output,
  input,
  mismatches = 3,
  threads = 1,
  parse_external = "-a -f",
  parse_internal = "a=TRUE, f=TRUE",
  import = "genome",
  keep_temp = FALSE,
  override = FALSE
)

Arguments

Details

Given a PAC object this function will extract the read sequences and align them against single or multiple fasta reference file(s). Summarized output will be saved as .Rdata files in the destination folder. In its default setting the function will only report hit or no hit in up to 3 mismatches, but this can easily be changed between 0-3 mismatches, to include alignment coordinates and reference name extractions. Information from the .Rdata files (Full_reanno_mis0/1/2/3.Rdata) can be extracted and added to a PAC object using make_reanno and add_reanno functions. For increased compatibility across platforms, the function provides both R internal Bowtie parsing (thorugh the Rbowtie package), as well as external parsing to a locally installed version of Bowtie.

Value

Will primarily generate .Rdata files in the destination folder (output) containing summarized information about the reference alignments. One file is generated for every mismatch specified in mismatches. The make_reanno function can then be used to extract and generate annotation tables for a PAC list object. Large temporary bowtie input and output files will also be generated in the destination folder, but are removed unless temp_remove=FALSE.

See Also

http://bowtie-bio.sourceforge.net/index.shtml for information about Bowtie and for Rbowtie: https://www.bioconductor.org/packages/release/bioc/html/Rbowtie.html. https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: PAC_annotate(), add_reanno(), as.reanno(), import_reanno(), make_conv(), make_reanno(), simplify_reanno()

Examples

######################################################### 
##### Simple example for reference mapping 
##### Please see manual for ?simply_reanno an ?add_reanno for more advanced 
##### mapping
closeAllConnections()
######################################################### 
##### Create an reanno object

### First, if you haven't already generated Bowtie indexes for the included
# fasta references you need to do so. If you are having problem see the small
# RNA guide (vignette) for more info.
                                                             
 ## tRNA:
 trna_file <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 trna_dir <- dirname(trna_file)
 
 if(!sum(stringr::str_count(list.files(trna_dir), ".ebwt")) ==6){     
     Rbowtie::bowtie_build(trna_file, 
                           outdir=trna_dir, 
                           prefix="tRNA", force=TRUE)
                           }
 ## rRNA:
 rrna_file <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 rrna_dir <- dirname(rrna_file)
 
 if(!sum(stringr::str_count(list.files(rrna_dir), ".ebwt")) ==6){
     Rbowtie::bowtie_build(rrna_file, 
                           outdir=rrna_dir, 
                           prefix="rRNA", force=TRUE)
                           }
                           
                           
##  Then load a PAC-object and remove previous mapping from anno:
 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 input <- list(trna= trna_file, rrna= rrna_file)

                                    
##  You may add an output path of your choice, but here we use a temp folder:
 output <- paste0(tempdir(),"/seqpac/test")


##  Then map the PAC-object against the fasta references. 
# Warning: if you use your own data, you may want to use override=FALSE, to avoid
# deleting previous mapping by mistake.

map_reanno(pac, input=input, output=output,
               type="internal", mismatches=0,  import="biotype", 
               threads=2, keep_temp=FALSE, override=TRUE)
 
##  Then import and generate a reanno-object of the temporary bowtie-files
reanno_biotype <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
 
                                    
## Now make some search terms against reference names to create shorter names
# Theses can be used to create factors in downstream analysis
# One search hit (regular expressions) gives one new short name 

bio_search <- list(
              rrna=c("5S", "5.8S", "12S", "16S", "18S", "28S", "pre_S45"),
              trna =c("_tRNA", "mt:tRNA"))
 
 
## You can merge directly with your PAC-object by adding your original 
# PAC-object, that you used with map_reanno, to merge_pac option.

pac <- add_reanno(reanno_biotype, bio_search=bio_search, 
                       type="biotype", bio_perfect=FALSE, 
                       mismatches = 0, merge_pac=pac)

table(anno(pac)$mis0_bio) 

############################################################################ 
## Similarly, you can use Bowtie indexed genome fasta references
## But don't forget to use import="genome" for coordinate import
#    
# input <- list(genome="<path_to_bowtie_indexed_fasta>")
# output_genome <- "<your_path_to_output_folder>"

## We can actually map against several genomes simultaneously:
# (to illustrate the principle we run the mycoplasma genome twice)

 mycoplasma_file <- system.file("extdata/mycoplasma_genome",
                                "mycoplasma.fa",
                                package = "seqpac", mustWork = TRUE)
 mycoplasma_dir<- gsub("mycoplasma.fa", "", mycoplasma_file)
 
## Don't forget to create bowtie indexes
 if(!sum(stringr::str_count(list.files(mycoplasma_dir), ".ebwt")) ==6){
     Rbowtie::bowtie_build(mycoplasma_file, 
                           outdir=mycoplasma_dir, 
                           prefix="mycoplasma", force=TRUE)
                           }
                           
 input <- list(genome1=mycoplasma_file, genome2=mycoplasma_file)
 map_reanno(PAC=pac, input=input, output=output, 
            type="internal", mismatches=0, import="genome", 
            threads=2, keep_temp=TRUE, override=TRUE)
 reanno_genome <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)
 
 table(overview(reanno_genome)$Any) # Low Mycoplasma contamination
 
############################################################################

Description

merge_lanes Identifies, reads and merges fastq lane files.

Usage

merge_lanes(input, output, threads = 1, nlanes = NULL)

Arguments

Details

Given an input path were flow cell lane files in fastq format (typically generated by for example Illumina sequencers) this function will try to identify unique samples among the lane files and merge all lane files for each sample. The function uses md5 checks to control that each lane file is unique (no accidental copies) and will compare expected outcomes given that all samples should have the same number of lanes (coming from the same experiment). If the automatic regex for commonly used lane naming does not work, try defining number of lanes to be merged by @param nlanes.

Value

Merged fastq files in destination folder.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), PAC_create(), make_PAC(), make_counts(), make_cutadapt(), make_pheno(), make_trim()

Examples

## The simple principle: 
# input <- "/some/path/to/lane/files/fastq.gz"
# output <- "/some/path/to/merged/files/"
# merge_lanes(input, output, threads=12)


## Real example
# First generate some correct file names (see above).
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
fq <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)
 
# Create an output folder
input <- paste0(tempdir(), "/lanes/")
output <- paste0(tempdir(), "/merged/")
dir.create(input, showWarnings=FALSE)
dir.create(output, showWarnings=FALSE)
# Fix compatible file names
file.copy(from = fq, to = input)
old_fls <- list.files(input, full.names=TRUE)
new_sample <- c(rep("sample1_", times=3), rep("sample2_", times=3))
new_lane <- rep(c("lane1","lane2","lane3"), times=2)
new_fls <- paste0(input,new_sample, new_lane, ".fastq.gz")
file.rename(from = old_fls, to = new_fls)

# Then merge the fastq files
merge_lanes(input, output, threads=2)

# You will find the files in:
input
output


##-----------------------------------------##
## Warning: Clean up temp folder           ##
# (Sometimes needed for automated examples) 

closeAllConnections()
fls_temp  <- tempdir()
fls_temp  <- list.files(fls_temp, recursive=TRUE, full.names = TRUE)
suppressWarnings(file.remove(fls_temp))

Description

names of objects in reanno

Usage

## S4 method for signature 'PAC'
names(x)

## S4 method for signature 'reanno'
names(x)

Arguments

Value

The names of the content in a PAC-object.

Names of the items in the reanno-object.

See Also

Other PAC methods: anno(), counts(), norm(), pheno(), show.PAC(), summary()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

######################################################### 
##### Create a reanno object

##  First load a PAC- object

 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 
##  Then specify paths to fasta references
# If you are having problem see the vignette small RNA guide for more info.
 
 trna_path <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)  
 rrna_path <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 
 ref_paths <- list(trna= trna_path, rrna= rrna_path)
                                    
##  Add output path of your choice.
# Here we use the R temporary folder depending on platform
if(grepl("windows", .Platform$OS.type)){
 output <- paste0(tempdir(), "\\seqpac\\test")
}else{
 output <- paste0(tempdir(), "/seqpac/test")}

## Make sure it is empty (otherwise you will be prompted for a question)
out_fls  <- list.files(output, recursive=TRUE)
closeAllConnections()
suppressWarnings(file.remove(paste(output, out_fls, sep="/")))

##  Then map your PAC-object against the fasta references
 map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=2,  import="biotype", 
               threads=2, keep_temp=FALSE)

##  Then generate a reanno-object of the temporary bowtie-files 
reanno_object <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)

## Accessing content and S4/S3 conversion:
names(reanno_object)
overview(reanno_object)
full(reanno_object)
rownames(reanno_object)
length(reanno_object)
nrow(reanno_object)
reanno_s3 <- as(reanno_object, "list")
reanno_s4 <- as.reanno(reanno_s3)

Description

number of samples

Usage

## S4 method for signature 'PAC'
ncol(x)

Arguments

Value

The number of samples in a PAC-object.

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

Access the norm list with normalized table in S4 PAC

Usage

norm(object)

## S4 method for signature 'PAC'
norm(object)

## S4 replacement method for signature 'PAC,list'
norm(object) <- value

Arguments

Value

The list of normalized tables (data.frames) from a PAC-object.

Updated S4 PAC-object.

See Also

Other PAC methods: anno(), counts(), names,PAC-method, pheno(), show.PAC(), summary()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

number of sequences

number of sequences

Usage

## S4 method for signature 'PAC'
nrow(x)

## S4 method for signature 'reanno'
nrow(x)

Arguments

Value

The number of sequences in a PAC-object.

Number of sequences in the reanno-object.

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

# extra fuctionality with s4 PAC-object:
names(pac)
length(pac)
nrow(pac)
ncol(pac)
rownames(pac)
colnames(pac)
pheno(pac) 
head(anno(pac))
head(counts(pac))
head(norm(pac)$cpm)

Description

Access the overview table in a S4 reanno object

Usage

overview(object)

## S4 method for signature 'reanno'
overview(object)

## S4 replacement method for signature 'reanno,data.frame'
overview(object) <- value

Arguments

Value

The overview table of a reanno-object as a tibble data.frame (class tbl_df, tbl, data.frame).

Updated reanno-object.

See Also

Other reanno methods: full()

Examples

######################################################### 
##### Create a reanno object

##  First load a PAC- object

 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))
 anno(pac) <- anno(pac)[,1, drop = FALSE]
 
 
##  Then specify paths to fasta references
# If you are having problem see the vignette small RNA guide for more info.
 
 trna_path <- system.file("extdata/trna", "tRNA.fa", 
                          package = "seqpac", mustWork = TRUE)  
 rrna_path <- system.file("extdata/rrna", "rRNA.fa", 
                          package = "seqpac", mustWork = TRUE)
 
 ref_paths <- list(trna= trna_path, rrna= rrna_path)
 
##  Add output path of your choice.
# Here we use the R temporary folder depending on platform 
if(grepl("windows", .Platform$OS.type)){
 output <- paste0(tempdir(), "\\seqpac\\test")
}else{
 output <- paste0(tempdir(), "/seqpac/test")}

## Make sure it is empty (otherwise you will be prompted for a question)
out_fls  <- list.files(output, recursive=TRUE)
closeAllConnections()
suppressWarnings(file.remove(paste(output, out_fls, sep="/")))

##  Then map your PAC-object against the fasta references 
 map_reanno(pac, input=ref_paths, output=output,
               type="internal", mismatches=2,  import="biotype", 
               threads=2, keep_temp=FALSE)

##  Then generate a reanno-object of the temporary bowtie-files 
reanno_object <- make_reanno(output, PAC=pac, mis_fasta_check = TRUE)

## Accessing content and S4/S3 conversion:
names(reanno_object)
overview(reanno_object)
full(reanno_object)
rownames(reanno_object)
length(reanno_object)
nrow(reanno_object)
reanno_s3 <- as(reanno_object, "list")
reanno_s4 <- as.reanno(reanno_s3)

Description

Given a PAC object and targets, this wrapper will perform several analyses of the PAC object and return a PDF with plots. More specifically, the plot will contain a jitterplot of fold changes, stack bar of mean counts/cpm, size distribution stack bars, piecharts, PCA plots, and if specified, a volcano plot based off adjusted p-values and fold changes from a DESeq2 analysis.

Usage

PAC_analyze(
  PAC,
  pheno_target = NULL,
  norm = NULL,
  anno_target = NULL,
  model = NULL
)

Arguments

Value

a list of plots in R, and a pdf file in home directory called "Results_seqpac.pdf".

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                 package = "seqpac", mustWork = TRUE))

                 
result_list <- PAC_analyze(pac, 
  pheno_target=list("stage"),
  norm="cpm", 
  anno_target=list("Biotypes_mis0"),
  model=~stage+batch)

Description

PAC_annotate Annotates the provided PAC object to a genome or biotypes

Usage

PAC_annotate(genome = NULL, biotype = NULL, output = NULL, PAC)

Arguments

Details

Given the path to either a reference genome, a biotype reference or both, PAC_annotate will wrap around all reanno-functions (map_reanno, make_reanno, and add_reanno) to annotate the PAC object. This wrapper comes with the caveat that it uses default values in most functions and is finetuned towards basic biotype annotation against the Ensembl ncRNA fasta reference. If that does not suit your analysis, we recommend you to run each function on its own for its full potential!

Value

PAC object with annotation information retrieved from mapping

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC reannotation: add_reanno(), as.reanno(), import_reanno(), make_conv(), make_reanno(), map_reanno(), simplify_reanno()

Examples

###########################################################

genome <- system.file("extdata/mycoplasma_genome", "mycoplasma.fa",
                     package = "seqpac", mustWork = TRUE)
genome_dir <- dirname(genome)

if(!sum(stringr::str_count(list.files(genome_dir), ".ebwt")) ==6){
  Rbowtie::bowtie_build(genome,
                        outdir=genome_dir,
                        prefix="mycoplasma", force=TRUE)
}

##  load a PAC-object and remove previous mapping from anno:
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata",
                 package = "seqpac", mustWork = TRUE))
anno(pac) <- anno(pac)[,1, drop = FALSE]

output <- paste0(tempdir(),"/seqpac/test")

#as we only test this against the mycoplasma genome,
#we expect a small overlap (0.03%)

pac <- PAC_annotate(PAC=pac,
                    genome = genome,
                    output = output)

Description

PAC_check Quickly checks rownames in PAC.

Usage

PAC_check(PAC, sample_check = NULL, seq_check = NULL)

Arguments

Details

Given a PAC object the function will cross check the rownames and colnames across Pheno, Anno and Counts.

Value

If PAC format is correct will return logical TRUE, else it will throw an error.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_create(), make_PAC(), make_counts(), make_cutadapt(), make_pheno(), make_trim(), merge_lanes()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

PAC_check(pac)  # TRUE

# Remove first row in Anno and check compatibility
new_Anno <- anno(pac)[-1,]
#PAC_check(pac, seq_check=new_Anno) # Error

# Add to pac an check
# anno(pac) <- new_Anno #error

Description

PAC_covplot Plotting sequences in a PAC object using an PAC mapping object.

Usage

PAC_covplot(
  PAC,
  map,
  summary_target = NULL,
  map_target = NULL,
  style = "line",
  xseq = TRUE,
  colors = NULL,
  check_overide = FALSE
)

Arguments

Details

Given a PAC object and a map object generated by PAC_mapper this function will attempt to plot the sequence coverage over long references.

Value

Plot list with plots generated by ggplot2

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

###########################################################
### Simple example of how to use PAC_mapper and PAC_covplot
# Note: More details, see vignette and manuals.)
# Also see: ?map_rangetype, ?tRNA_class or ?PAC_trna for more examples
# on how to use PAC_mapper.

## Load PAC-object, make summaries and extract rRNA and tRNA
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))

pac <- PAC_summary(pac, norm = "cpm", type = "means", 
                   pheno_target=list("stage", unique(pheno(pac)$stage)))

pac_tRNA <- PAC_filter(pac, anno_target = list("Biotypes_mis0", "tRNA"))


## Give paths to a fasta reference (with or without bowtie index)
#  (Here we use an rRNA/tRNA fasta included in seqpac) 

ref_tRNA <- system.file("extdata/trna", "tRNA.fa", 
                         package = "seqpac", mustWork = TRUE)


## Map using PAC-mapper
      
map_tRNA <- PAC_mapper(pac_tRNA, mismatches=0, 
                        threads=1, input=ref_tRNA, override=TRUE)

## Plot tRNA using xseq=TRUE gives you reference sequence as X-axis:
# (OBS! Long reference will not print well with xseq=TRUE)
cov_tRNA <- PAC_covplot(pac_tRNA, map_tRNA, 
                        summary_target = list("cpmMeans_stage"), 
                        xseq=TRUE, style="line", 
                        color=c("red", "black", "blue"))

cov_tRNA[[1]]

## Check which tRNA reached a decent number of unique fragments 
# (OBS! This is a very down sampled dataset)
logi_hi <- unlist(lapply(map_tRNA, function(x){nrow(x$Alignments) > 10 }))

table(logi_hi)  
names(map_tRNA[logi_hi])

targets <- c("Ala-AGC-1-1", "Lys-CTT-1-13","Ser-AGA-2-2")

cov_tRNA_sub <- PAC_covplot(pac_tRNA, map_tRNA, 
                        summary_target = list("cpmMeans_stage"),
                        map_target = targets,
                        xseq=TRUE, style="line", 
                        color=c("red", "black", "blue"))

cowplot::plot_grid(plotlist= cov_tRNA_sub)

Description

PAC_create Merges lanes, performs adaptor trimming and makes a first PAC object with count and basic phenotypic information

Usage

PAC_create(lanes = NULL, trim = NULL, input, output = NULL)

Arguments

Details

Given (at minimum) an path to input files to be analyzed, PAC_create will wrap around all necessary preparatory steps (merge_lanes, make_counts, and make_PAC) to produce a PAC object. This wrapper comes with the caveat that it uses default values in most functions. If that does not suit your analysis, we recommend you to run each function on its own!

Value

PAC object with values in pheno and count

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC generation: PAC_check(), make_PAC(), make_counts(), make_cutadapt(), make_pheno(), make_trim(), merge_lanes()

Examples

###########################################################
### test the map_rangetype function
# More complicated examples can be found in the vignette.
##----------------------------------------

# First create an annotation blank PAC with group means
sys_path = system.file("extdata", package = "seqpac", mustWork = TRUE)
fq <- list.files(path = sys_path, pattern = "fastq", all.files = FALSE,
                full.names = TRUE)
 
# Create an output folder
input <- paste0(tempdir(), "/lanes/")
output <- paste0(tempdir(), "/merged/")
dir.create(input, showWarnings=FALSE)
dir.create(output, showWarnings=FALSE)

# Fix compatible file names
file.copy(from = fq, to = input)
old_fls <- list.files(input, full.names=TRUE)
new_sample <- c(rep("sample1_", times=3), rep("sample2_", times=3))
new_lane <- rep(c("lane1","lane2","lane3"), times=2)
new_fls <- paste0(input,new_sample, new_lane, ".fastq.gz")
file.rename(from = old_fls, to = new_fls)

# Then merge the fastq files
pac <- PAC_create(lanes=TRUE, trim="default_neb", input, output)

Description

PAC_deseq DESeq2 analysis on PAC_object.

Usage

PAC_deseq(
  PAC,
  model,
  deseq_norm = FALSE,
  test = "Wald",
  fitType = "local",
  threads = 1,
  pheno_target = NULL
)

Arguments

Details

Given a PAC object this function will apply a differential expression analysis using DESeq2.

Value

A list of objects:

1st object - Summarized result table merged with anno(PAC)

2nd object - Target graphs (p-val distribution and volcano)

3rd object - All output from DESeq

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

# Note, these examples will generate some warnings since data is based on
# heavily down-sampled fastq files, where many sequences receives low counts
# in specific groups.

## Load test data
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                package = "seqpac", mustWork = TRUE))

## Simple model with embryonic stages using Wald test with local fit (default)
table(pheno(pac)$stage)
output_deseq <- suppressWarnings(PAC_deseq(pac, model= ~stage, threads=1))

## Batch corrected, graphs are generated for 'stage' (=first in the model)
output_deseq <- suppressWarnings(PAC_deseq(pac, model= ~stage + batch,
                                          threads=1))

## Using pheno_target 
output_deseq <- suppressWarnings(PAC_deseq(pac,model= ~stage + batch, 
                                          pheno_target=list("batch"),
                                          threads=1))

## With pheno_target we can change the direction for the comparison
# Stage5 vs Stage3 (reverse order):
output_deseq <- suppressWarnings(PAC_deseq(pac, model= ~stage + batch, 
                         pheno_target = list("stage", c("Stage5", "Stage3")),
                         threads=1))  

## In the output you find PAC merged results, target plots and output_deseq   
names(output_deseq)
tibble::as_tibble(output_deseq$result)

Description

PAC_filter Filter PAC objects.

Usage

PAC_filter(
  PAC,
  nucleotide_range = NULL,
  threshold = 0,
  coverage = 0,
  norm = "counts",
  subset_only = FALSE,
  stat = FALSE,
  pheno_target = NULL,
  anno_target = NULL,
  allbut = FALSE
)

Arguments

Details

Given a PAC object the function will extract sequences within a given size interval and percent coverage across independent samples.

Value

A list of objects: PAC object with filtered data. (optional) A coverage plot

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                 package = "seqpac", mustWork = TRUE))

###--------------------------------------------------------------------- 
## Extracts all sequences between 10-80 nt in length with at 
## least 5 counts in 20% of all samples.
pac_lowfilt <- PAC_filter(pac, nucleotide_range=c(10,80), threshold=5,
                         coverage=20, norm = "counts",
                         pheno_target=NULL, anno_target=NULL)

###--------------------------------------------------------------------- 
## Extracts sequences with 22 nt size and the samples in Batch1 and Batch2.
pac_subset <- PAC_filter(pac, subset_only = TRUE,
                        pheno_target=list("batch", c("Batch1", "Batch2")), 
                        anno_target=list("Size", "22"))

###--------------------------------------------------------------------- 
## Extracts all sequences with >=5 counts in 100% of samples a within stage
filtsep <- PAC_filtsep(pac, norm="counts", threshold=5, 
                      coverage=100, pheno_target= list("stage"))

pac_filt <- PAC_filter(pac, subset_only = TRUE,
                     anno_target= unique(do.call("c", as.list(filtsep))))
                     
###--------------------------------------------------------------------- 
## Extract all biotypes but tRNA in all samples but in sample fastq_1

pac_filt <- PAC_filter(pac, anno_target=list("Biotypes_mis0","tRNA"),
                     pheno_target=list("sample","fastq1"), allbut=TRUE)

Description

PAC_filtsep Group separated filtered name list

Usage

PAC_filtsep(
  PAC,
  norm = "counts",
  threshold = 10,
  coverage = 100,
  pheno_target = NULL,
  output = "sequence"
)

Arguments

Details

Given a PAC object the function will extract sequences within the given filter within a given group.

Value

A data.frame (see output for details).

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_gtf(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))

## Keep sequences with 5 counts (threshold) in 100% (coverage) of 
## samples in a group:
 # Use PAC_filtsep to find sequences 
 filtsep <- PAC_filtsep(pac, norm="counts", threshold=5, 
                        coverage=100, pheno_target= list("stage"))
                        
 # Filter by unique sequences passing filtsep  
 filtsep <- unique(do.call("c", as.list(filtsep)))
 pac_filt <- PAC_filter(pac, subset_only = TRUE, anno_target= filtsep)
 
 # Find overlap
 olap <- reshape2::melt(filtsep, 
                        measure.vars = c("Stage1", "Stage3", "Stage5"), 
                        na.rm=TRUE)
                        
## Upset plot using the UpSetR package
 # (when output="binary" PAC_filtsep uses filtsep_bin for binary conversion
 # Use PAC_filtsep with binary output
 filtsep_bin <- PAC_filtsep(pac, norm="counts", threshold=5, 
                            coverage=100, pheno_target= list("stage"), 
                            output="binary")
 
# Plot Wenn diagram or UpSetR
#
# plot(venneuler::venneuler(data.frame(olap[,2], olap[,1]))) 
#
# UpSetR::upset(filtsep_bin, sets = colnames(filtsep_bin), 
#              mb.ratio = c(0.55, 0.45), order.by = "freq", keep.order=TRUE)

Description

This function will annotate a PAC object using input from a GTF/GFF file.

Usage

PAC_gtf(
  PAC,
  genome = NULL,
  mismatches = 3,
  return = "simplify",
  stranded = FALSE,
  gtf = NULL,
  targets = NULL,
  threads = 1
)

Arguments

Details

Given a PAC object and a gtf formatted annotation file(s), this function will attempt to annotate sequences mapped to a reference genome against genomic coordinates in the gtf file(s). In case no genomic mapping coordinates are available in the PAC object, the function provides a backdoor into the PAC reannotation workflow, where genome mapping is performed using Bowtie.

Value

List, tibble dataframe or updated PAC object. See option return for more information.

See Also

https://github.com/OestLab/seqpac for updates.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_jitter(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

# Load PAC
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                 package = "seqpac", mustWork = TRUE))

# Create a gtf file from PAC coordinates
anno <- anno(pac)
anno <- anno[!grepl("Warning", anno$mis0_chromosomes_genome),]
anno <- anno[!is.na(anno$mis0_chromosomes_genome),]
coord <- anno$mis0_chromosomes_genome
coord <- suppressWarnings(do.call("rbind", strsplit(coord, "\\|"))[,1])
coord <- suppressWarnings(do.call("rbind", strsplit(coord, "\\;start=|;")))
gr <- GenomicRanges::GRanges(seqnames=coord[,1], 
                             IRanges::IRanges(as.numeric(coord[,2]), 
                                              as.numeric(coord[,2])+
                                              anno$Size ), 
                             strand=coord[,3])

GenomicRanges::mcols(gr) <- data.frame(biotype=anno$Biotypes_mis3, 
                                       bio_zero=as.character(anno$mis0_bio))
spl <- sample(1:length(gr), round(length(gr)*0.3), replace=FALSE)
gr1 <- gr[spl]
gr2 <- gr[!1:length(gr) %in% spl]

# Prepare temp folder and save artifical gtf
out1 <- paste0(tempdir(), "/temp1.gtf")
out2 <- paste0(tempdir(), "/temp2.gtf")

rtracklayer::export(gr1, out1, format="gtf")
rtracklayer::export(gr2, out2, format="gtf")

# Make sure PAC contains full genome coordinates
# (In the add_reanno function you may choose how many coordinates to report)
# (If there are more, a 'Warning' will be added to the annotation)
# (Here we remove those to avoid problems)

new_anno <-  anno(pac) [, grepl("chromosomes_genome|Size", 
                               colnames(anno(pac)))]
test <- new_anno[, 3:ncol(new_anno)]
test <- apply(test, 2, function(x){substr(x, 1, 7)})
test <- apply(test, 1, function(x){paste(x, collapse = "")})
new_anno$temp <- ifelse(grepl("Warning",  test), "rm", "keep")
anno(pac) <- new_anno
pac_small <- PAC_filter(pac, subset_only=TRUE, nucleotide_range = c(20,22),
                  anno_target=list("temp", "keep"))

#  Run PAC_gtf
gtf <- list(gtf1=out1, gtf2=out2)
target <- list(gtf1=c("biotype","bio_zero"), gtf2=c("biotype","bio_zero"))

pac_merge <- PAC_gtf(pac_small, mismatches=0, return="merge", 
                    gtf=gtf, target=target, threads=2)

Description

PAC_jitter Plots a jitter plot using the information in the Anno and Counts tables in a PAC object.

Usage

PAC_jitter(
  PAC,
  summary_target = NULL,
  anno_target = NULL,
  style = "jitter",
  limits = NULL,
  ypos_n = NULL,
  colors = NULL,
  box = TRUE
)

Arguments

Details

Given a PAC object with grouped summaries, the function will use column(s) in the Anno object to group assigned summary columns (summary_target) and plot a jitter plot.

Value

A plot-list object with jitter plots generated by ggplot2.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_mapper(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

## Prepare
load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                  package = "seqpac", mustWork = TRUE))
pac <- PAC_norm(pac, norm="cpm")

pac <- PAC_summary(PAC=pac, norm = "cpm", type = "log2FC", 
                   pheno_target=list("stage"))
pac <- PAC_summary(PAC=pac, norm = "cpm", type = "percentgrand", 
                   pheno_target=list("stage"))

## Jitter plots
plots_FC <- PAC_jitter(pac, summary_target=list("Log2FC_stage"), 
                       anno_target=list("Biotypes_mis0"))
plots_FCgrand <- PAC_jitter(pac, summary_target=list("percGrand_stage"), 
                            anno_target=list("Biotypes_mis0"))

cowplot::plot_grid(plotlist=plots_FC, nrow = 3, ncol = 1)
cowplot::plot_grid(plotlist=plots_FCgrand, nrow = 3, ncol = 1) 

## Violin plots instead
plots_FC <- PAC_jitter(pac, style="violin", 
                       summary_target=list("Log2FC_stage"), 
                       anno_target=list("Biotypes_mis0"))
plots_FCgrand <- PAC_jitter(pac, style="violin", 
                            summary_target=list("percGrand_stage"), 
                            anno_target=list("Biotypes_mis0"))

cowplot::plot_grid(plotlist=plots_FC, nrow = 3, ncol = 1)
cowplot::plot_grid(plotlist=plots_FCgrand, nrow = 3, ncol = 1)  


## Violin with changed biotype order
new_order  <- as.character(unique(anno(pac)$Biotypes_mis0))[c(2,4,3,6,7,5,1)]
plots_FC <- PAC_jitter(pac, style="violin", 
                       summary_target=list("Log2FC_stage"), 
                       anno_target=list("Biotypes_mis0", new_order))
cowplot::plot_grid(plotlist=plots_FC, nrow = 3, ncol = 1)

Description

PAC_mapper Map sequences against a small reference.

Usage

PAC_mapper(
  PAC,
  input,
  mismatches = 0,
  multi = "remove",
  threads = 1,
  N_up = "",
  N_down = "",
  report_string = FALSE,
  override = FALSE
)

Arguments

Details

Given a PAC object and the path to a fasta reference file, this function will map sequences in the PAC using a 'backdoor' into the reanno workflow.

Value

Stacked list, where each object on the highest level contains: (Object 1) Reference name and sequence. (Object 2) Data.frame showing the mapping results of each query sequence that mapped to Object 1.

See Also

https://github.com/OestLab/seqpac for updates on the current package.

Other PAC analysis: PAC_covplot(), PAC_deseq(), PAC_filter(), PAC_filtsep(), PAC_gtf(), PAC_jitter(), PAC_nbias(), PAC_norm(), PAC_pca(), PAC_pie(), PAC_saturation(), PAC_sizedist(), PAC_stackbar(), PAC_summary(), PAC_trna(), as.PAC(), filtsep_bin(), map_rangetype(), tRNA_class()

Examples

###########################################################
### Simple example of how to use PAC_mapper 
# Note: More details, see vignette and manuals.)
# Also see: ?map_rangetype, ?tRNA_class or ?PAC_trna, ?PAC_covplot 
# for more examples on how to use PAC_mapper.

## Load PAC-object, make summaries and extract rRNA and tRNA
 load(system.file("extdata", "drosophila_sRNA_pac_filt_anno.Rdata", 
                   package = "seqpac", mustWork = TRUE))

pac <- PAC_summary(pac, norm = "cpm", type = "means", 
                   pheno_target=list("stage", unique(pheno(pac)$stage)))
                   
pac_rRNA <- PAC_filter(pac, anno_target = list("Biotypes_mis0", "rRNA"))

## Give paths to a fasta reference (with or without bowtie index)
#  (Here we use an rRNA/tRNA fasta included in seqpac) 

ref_rRNA <- system.file("extdata/rrna", "rRNA.fa", 
                         package = "seqpac", mustWork = TRUE)
                         

## Map using PAC-mapper
map_rRNA <- PAC_mapper(pac_rRNA, mismatches=0, 
                        threads=1, input=ref_rRNA, override=TRUE)

Description

PAC_nbias analyses nucleotide bias.

Usage

PAC_nbias(
  PAC,
  position = 1,
  norm = NULL,
  nucleotide_range = NULL,
  anno_target = NULL,
  pheno_target = NULL,
  summary_target = NULL,
  colors = NULL,
  ymax = NULL,
  data_only = FALSE
)

Arguments