Pipeline for Bulk Data

Semi-supervised isofrom detection and annotation for long read data. This variant is meant for bulk samples. Specific parameters relating to analysis can be changed either through function arguments, or through a configuration JSON file.

Usage

bulk_long_pipeline(
  annotation,
  fastq,
  outdir,
  genome_fa,
  minimap2 = NULL,
  k8 = NULL,
  config_file = NULL
)

Arguments

annotation: The file path to the annotation file in GFF3 format
fastq: The file path to input fastq file
outdir: The path to directory to store all output files.
genome_fa: The file path to genome fasta file.
minimap2: Path to minimap2, if it is not in PATH. Only required if either or both of do_genome_align and do_read_realign are TRUE.
k8: Path to the k8 Javascript shell binary. Only required if do_genome_align is TRUE.
config_file: File path to the JSON configuration file. If specified, config_file overrides all configuration parameters

Value

if do_transcript_quantification set to true, bulk_long_pipeline returns a SummarizedExperiment object, containing a count matrix as an assay, gene annotations under metadata, as well as a list of the other output files generated by the pipeline. The pipeline also outputs a number of output files into the given outdir directory. These output files generated by the pipeline are:

transcript_count.csv.gz: - a transcript count matrix (also contained in the SummarizedExperiment)
isoform_annotated.filtered.gff3: - isoforms in gff3 format (also contained in the SummarizedExperiment)
transcript_assembly.fa: - transcript sequence from the isoforms
align2genome.bam: - sorted BAM file with reads aligned to genome
realign2transcript.bam: - sorted realigned BAM file using the transcript_assembly.fa as reference
tss_tes.bedgraph: - TSS TES enrichment for all reads (for QC)

if do_transcript_quantification set to false, nothing will be returned

Details

By default FLAMES use minimap2 for read alignment. After the genome alignment step (do_genome_align), FLAMES summarizes the alignment for each read by grouping reads with similar splice junctions to get a raw isoform annotation (do_isoform_id). The raw isoform annotation is compared against the reference annotation to correct potential splice site and transcript start/end errors. Transcripts that have similar splice junctions and transcript start/end to the reference transcript are merged with the reference. This process will also collapse isoforms that are likely to be truncated transcripts. If isoform_id_bambu is set to TRUE, bambu::bambu will be used to generate the updated annotations. Next is the read realignment step (do_read_realign), where the sequence of each transcript from the update annotation is extracted, and the reads are realigned to this updated transcript_assembly.fa by minimap2. The transcripts with only a few full-length aligned reads are discarded. The reads are assigned to transcripts based on both alignment score, fractions of reads aligned and transcript coverage. Reads that cannot be uniquely assigned to transcripts or have low transcript coverage are discarded. The UMI transcript count matrix is generated by collapsing the reads with the same UMI in a similar way to what is done for short-read scRNA-seq data, but allowing for an edit distance of up to 2 by default. Most of the parameters, such as the minimal distance to splice site and minimal percentage of transcript coverage can be modified by the JSON configuration file (config_file).

The default parameters can be changed either through the function arguments are through the configuration JSON file config_file. the pipeline_parameters section specifies which steps are to be executed in the pipeline - by default, all steps are executed. The isoform_parameters section affects isoform detection - key parameters include:

Min_sup_cnt: which causes transcripts with less reads aligned than it's value to be discarded
MAX_TS_DIST: which merges transcripts with the same intron chain and TSS/TES distace less than MAX_TS_DIST
strand_specific: which specifies if reads are in the same strand as the mRNA (1), or the reverse complemented (-1) or not strand specific (0), which results in strand information being based on reference annotation.

Examples

# download the two fastq files, move them to a folder to be merged together
temp_path <- tempfile()
bfc <- BiocFileCache::BiocFileCache(temp_path, ask = FALSE)
file_url <-
  "https://raw.githubusercontent.com/OliverVoogd/FLAMESData/master/data"
# download the required fastq files, and move them to new folder
fastq1 <- bfc[[names(BiocFileCache::bfcadd(bfc, "Fastq1", paste(file_url, "fastq/sample1.fastq.gz", sep = "/")))]]
fastq2 <- bfc[[names(BiocFileCache::bfcadd(bfc, "Fastq2", paste(file_url, "fastq/sample2.fastq.gz", sep = "/")))]]
annotation <- bfc[[names(BiocFileCache::bfcadd(bfc, "annot.gtf", paste(file_url, "SIRV_isoforms_multi-fasta-annotation_C_170612a.gtf", sep = "/")))]]
genome_fa <- bfc[[names(BiocFileCache::bfcadd(bfc, "genome.fa", paste(file_url, "SIRV_isoforms_multi-fasta_170612a.fasta", sep = "/")))]]
fastq_dir <- paste(temp_path, "fastq_dir", sep = "/") # the downloaded fastq files need to be in a directory to be merged together
dir.create(fastq_dir)
file.copy(c(fastq1, fastq2), fastq_dir)
#> [1] TRUE TRUE
unlink(c(fastq1, fastq2)) # the original files can be deleted

outdir <- tempfile()
dir.create(outdir)
se <- bulk_long_pipeline(
  annotation = annotation, fastq = fastq_dir, outdir = outdir, genome_fa = genome_fa,
  config_file = create_config(outdir, type = "sc_3end", threads = 1, no_flank = TRUE)
)
#> Writing configuration parameters to:  /tmp/RtmpoS7Kzz/file1f6d2cd15227/config_file_8045.json 
#> #### Input parameters:
#> {
#>   "pipeline_parameters": {
#>     "seed": [2022],
#>     "threads": [1],
#>     "do_barcode_demultiplex": [true],
#>     "do_gene_quantification": [true],
#>     "do_genome_alignment": [true],
#>     "do_isoform_identification": [true],
#>     "bambu_isoform_identification": [false],
#>     "multithread_isoform_identification": [false],
#>     "do_read_realignment": [true],
#>     "do_transcript_quantification": [true],
#>     "oarfish_quantification": [true]
#>   },
#>   "barcode_parameters": {
#>     "max_bc_editdistance": [2],
#>     "max_flank_editdistance": [8],
#>     "pattern": {
#>       "primer": ["CTACACGACGCTCTTCCGATCT"],
#>       "BC": ["NNNNNNNNNNNNNNNN"],
#>       "UMI": ["NNNNNNNNNNNN"],
#>       "polyT": ["TTTTTTTTT"]
#>     },
#>     "strand": ["-"],
#>     "TSO_seq": ["AAGCAGTGGTATCAACGCAGAGTACATGGG"],
#>     "TSO_prime": [3],
#>     "full_length_only": [false]
#>   },
#>   "isoform_parameters": {
#>     "generate_raw_isoform": [false],
#>     "max_dist": [10],
#>     "max_ts_dist": [100],
#>     "max_splice_match_dist": [10],
#>     "min_fl_exon_len": [40],
#>     "max_site_per_splice": [3],
#>     "min_sup_cnt": [5],
#>     "min_cnt_pct": [0.001],
#>     "min_sup_pct": [0.2],
#>     "bambu_trust_reference": [true],
#>     "strand_specific": [0],
#>     "remove_incomp_reads": [4],
#>     "downsample_ratio": [1]
#>   },
#>   "alignment_parameters": {
#>     "use_junctions": [true],
#>     "no_flank": [true]
#>   },
#>   "realign_parameters": {
#>     "use_annotation": [true]
#>   },
#>   "transcript_counting": {
#>     "min_tr_coverage": [0.4],
#>     "min_read_coverage": [0.4]
#>   }
#> } 
#> gene annotation: /tmp/RtmpoS7Kzz/file1f6d5b8aad1/1f6d27d79d75_SIRV_isoforms_multi-fasta-annotation_C_170612a.gtf 
#> genome fasta: /tmp/RtmpoS7Kzz/file1f6d5b8aad1/1f6d60187dc6_SIRV_isoforms_multi-fasta_170612a.fasta 
#> input fastq files: /tmp/RtmpoS7Kzz/file1f6d5b8aad1/fastq_dir/1f6d23439419_sample1.fastq.gz
#>  /tmp/RtmpoS7Kzz/file1f6d5b8aad1/fastq_dir/1f6d31effc3d_sample2.fastq.gz
#> output directory: /tmp/RtmpoS7Kzz/file1f6d2cd15227 
#> minimap2 path: 
#> k8 path: 
#> #### Aligning reads to genome using minimap2
#> 	Aligning sample  1f6d23439419_sample1 ...
#> 02:12:50 AM Fri Oct 25 2024 minimap2_align
#> 	Aligning sample  1f6d31effc3d_sample2 ...
#> 02:12:53 AM Fri Oct 25 2024 minimap2_align
#> 02:12:56 AM Fri Oct 25 2024 find_isoform
#> #### Realign to transcript using minimap2
#> 	Realigning sample  1f6d23439419_sample1 ...
#> 02:12:57 AM Fri Oct 25 2024 minimap2_realign
#> file renamed to  1f6d23439419_sample1_realign2transcript.bam 
#> Warning: cannot remove file '/tmp/RtmpoS7Kzz/file1f6d2cd15227/1f6d23439419_sample1_tmp_align.bam', reason 'No such file or directory'
#> 	Realigning sample  1f6d31effc3d_sample2 ...
#> 02:12:58 AM Fri Oct 25 2024 minimap2_realign
#> file renamed to  1f6d31effc3d_sample2_realign2transcript.bam 
#> Warning: cannot remove file '/tmp/RtmpoS7Kzz/file1f6d2cd15227/1f6d31effc3d_sample2_tmp_align.bam', reason 'No such file or directory'
#> #### Generating transcript count matrix

Usage

Arguments

Value

Details

See also

Examples