1. Preliminary analysis - Bulk

1. Preliminary analysis - Bulk

The pre-processing workflows extract clean s-reads using 3 tools as below which can then be provided to the alignment applications and other downstream workflows similar to those used to analyze Isoseq data.

• skera for de-concatenating the MAS arrays into individual cDNA molecules and generate segmented reads (s-reads),

• lima to de-multiplex a bulk pool and remove unwanted combinations while orienting sequences in 5’ → 3’ orientation and

• isoseq refine for trimming poly(A) tails and extracting Full length non-concatemer reads (FLNC) from s-reads.

• samtools merge - only applicable in case technical replicates are generated to boost yeild by sequencing over mulitple flowcells, in which case the reads for each biological sample would require to me merged together.

1.1. pbskera

The pbskera workflow, as detailed below, processes raw HiFi reads generated with Sequel2e and Revio Long Read sequencers. The HiFi reads are a current default, and can be plugged in directly into the workflow to get segmented s-reads.

Workflow configuration for runnning these over cloud platforms supporting Cromwell like Terra can be found here:-

Dockstore : skera_w_QCplots.wdl

Github: Kinnex Preliminary Processing

Test Data can be found here (public, requester-pays) : gs://mdl-preprocess-refs/test_data/m12345_123456_123456_s1.hifi_reads.bcM0003.bam

The direct command executed here is:

    skera split –j 8 reads.hifi.bam mas_adapters.fasta reads.skera.bam

Input arguments for pbskera_main

Table 1 bulk - skera

Option name	example value	description
input_bam	this.input_bam	Replace “input_bam” with the column that contains the path to the input file in the sample table.
arraysize	8	Size of the MAS array, could be 8, 16 or 12 depending upon the library type
mas_adapters_fasta	“gs://mdl-preprocess-refs/MAS_adapters/mas8_primers.fasta”	MAS adapters FASTA file for de-concatenation. Ref files can be found at gs://mdl-preprocess-refs/MAS_adapters/
sample_id	this.movie_name	Replace “sample_id” with the column that contains the name of the samples in the sample table. Recommended id is the movie name
gcs_output_dir	“gs://output_dir/”	Output directory to organize intermediates and QC plots
num_threads	8	Number of threads (set atleast equal to the number of cpu)
cpu	8	Number of cpus

Running the workflow above also generates plots to sanity check the s-reads generated by deconacatenation. There are 3 plots as of now:

Fig. 1 Read Length Histogram

Fig. 2 Ligation Heatmap

Fig. 3 Concatenation Histogram

For a bulk 8-mer we expect the plots to be similar to as above, with maximum number of reads assigned to a complete 8-mer configuration. In addition, to the readlength plot, the concatenation histogram should also indicate high percentages (>90%) to be assigned to a concatenation factor of 8. The ligation heatmap distributes the number of reads by adapter pairs found in the array. They should cleanly align along the diagonal for a well-performing array.

1.2. lima demux + isoseq refine

This workflow uses 2 tools to extract clean s-reads from the skera.bam received above.

Workflow configuration for runnning these over cloud platforms supporting Cromwell like Terra can be found here:-

Dockstore : bulk_demux.wdl

Github : Kinnex Preliminary Processing

Test Data: gs://mdl-preprocess-refs/test_data/m12345_123456_123456_s1.hifi_reads.bcM0003.bam (public, requester-pays)

The direct command executed here is:

    lima --isoseq --log-level INFO -j 8 reads.skera.bam bulk_barcodes.fasta sample_id.lima.bam

    isoseq refine --require-polya -j 8 sample_id.lima.bam bulk_barcodes.fasta sample_id.refine.bam

Input arguments for bulk_demux

Table 2 bulk lima+refine

Option name	example value	description
input_bam	this.skera_bam	Replace “skera_bam” with the datatable column or a gs://path to s-reads generated with skera in part1.
sample_id	this.sample_id	Replace “sample_id” with the column or the sample name.
trimPolyA	“True”	Boolean set as “True”/”False” to remove polyA tails
clipAdapters	“True”	Boolean set as “True”/”False” to remove MAS adapters
bulk_barcodes_fasta	“gs://mdl-preprocess-refs/Bulk_Lima_barcodes/bulk12_barcodes.fasta”	Provide bulk_barcodes.fasta corresponding to the single cell kit used in library prep. Ref files can be found at
gcs_output_dir	“gs://output_dir/”	Output directory to organize intermediates and QC plots
num_threads	16	Number of threads (number of cpus are set equal to the number of threads)
Outputs:
lima_out	“gs://output_dir/lima”	Output directory for lima outs
refine_out	“gs://output_dir/refine”	Output directory for refine outs

1.3. merge

The merge_refine_bams workflow, as detailed below, merges refined clean s-reads generated with workflows above for various technical replicates. It is an optional workflow and can be replaced by merging replicates by hand which is recommended if custom edits are to be made before merging like adding custom tags to the bam for each replicate.

Workflow configuration for runnning these over cloud platforms supporting Cromwell like Terra can be found here:-

Dockstore : merge_refine_bams.wdl

Github: Kinnex Preliminary Processing

Test Data can be found here (public, requester-pays) : gs://mdl-preprocess-refs/test_data/m12345_123456_123456_s1.hifi_reads.bcM0003.bam

It requires an inputfile.csv with the sample to barcode combination used for multiplexing the samples. An example file corresponding to the test data is provided below.`Sample_ID`,`Kinnex_Adapter`,`IsoSeq_primer` are required (case-insensitive) columns. The file can have additional optional columns with snake_cased colnames.

Table 3 bulk - idmap.csv , sample to barcode matching

Sample_ID	Kinnex_Adapter	IsoSeq_primer
test_bcM0001_bc01	bcM0001	bc01
test_bcM0001_bc02	bcM0001	bc02
test_bcM0001_bc03	bcM0001	bc03
test_bcM0001_bc04	bcM0001	bc04
test_bcM0001_bc05	bcM0001	bc05
test_bcM0001_bc06	bcM0001	bc06
test_bcM0001_bc07	bcM0001	bc07
test_bcM0001_bc08	bcM0001	bc08
test_bcM0001_bc09	bcM0001	bc09
test_bcM0001_bc10	bcM0001	bc10
test_bcM0001_bc11	bcM0001	bc11
test_bcM0001_bc12	bcM0001	bc12

The workflow calls the script below which leverages pysam utils to merge reads and rename by the sample_ids provided in the idmap.csv file above the here is:

   python /usr/local/src/masseq_data_processing/pb_plots/mergeBams.py \
         -idmap ~{barcode_to_sample} \
         -bampath . \
         -limacountsdir . \
         -outdir . \
         -mergeReplicates \
         -setTitleSamplePlot ~{datasetId}

Input arguments for merge_refine_bams

Table 4 bulk - merge

option	example value	description
barcode_to_sample	“gs://kinnex_example_data/bulk8_test_data/ kinnex_bulk8_test_idmap.csv”	csv with sampleid to barcode mapping as described above
bulk_barcodes_fasta	“gs://mdl-preprocess-refs/Bulk_Lima_barcodes/ bulk12_barcodes.fasta”
datasetId	“kinnex_bulk8_test”
gcs_output_dir	“gs://bucket/out_path”	GS location where merge folder will be created to hold merged reads and QC tables and plot
lima_out	“gs://kinnex_example_data/bulk8_test_data/lima”
refine_out	“gs://kinnex_example_data/bulk8_test_data/refine”
mergeBams	“True”	Boolean value default “True” to merge, “False” for only generate counts plot
cpu	“${8}”
num_threads	“${8}”

Upon merging the replicates, the reads are ready to be analysed with the Secondary Processing workflows.