Usage questions

[olgabot]

Hello Liz,
Hope you are doing well :) I'm using your latest version (b7ca1f3) and I had some questions. I'm quite new to gene therapy/rAAVs so forgive me if these questions are rudimentary.

1. What reads "should be coming from the reference genome"?

In the tutorial, I'm confused about what "expected range of where the reads should be coming from in the reference genome" means, and how that is different from the ITR-to-ITR insert specified in the annotation.txt file? Is this where random bits of DNA could get incorporated into the 5' ITR -- gene of interest -- 3' ITR product?

For the ssAAV aligned BAM file, we can specify an expected range of where the reads should be coming from in the reference genome
(emphasis mine)

2. plotAAVreport.R doesn't seem up to date

Also, it seems the plotAAVreport.R script has not yet been updated since the changes made on 6/30, because when I used the prefix as specified in the tutorial, the script seems to be opening the folder, rather than looking for files with that prefix:

`Lot-PB19045--summarize-aav-alignment` files

(gene-therapy)
 Tue  5 Jul - 23:38  /data/bbio-gt/2022-06-30--raav-seq-pilot/CCS 
 olgabot@x86_64-conda-linux-gnu  ll Lot-PB19045--summarize-aav-alignmen*
Permissions Size User    Date Modified Name
.rw-rw-r--   59M olgabot  5 Jul 19:09  Lot-PB19045--summarize-aav-alignment.nonmatch_stat.csv
.rw-rw-r--  4.0M olgabot  5 Jul 19:09  Lot-PB19045--summarize-aav-alignment.per_read.csv
.rw-rw-r--   11M olgabot  5 Jul 19:09  Lot-PB19045--summarize-aav-alignment.summary.csv
.rw-rw-r--   57M olgabot  5 Jul 19:09  Lot-PB19045--summarize-aav-alignment.tagged.bam

Here's the error message:

(gene-therapy)
 Tue  5 Jul - 19:22  /data/bbio-gt/2022-06-30--raav-seq-pilot/CCS 
 olgabot@x86_64-conda-linux-gnu  Rscript ~/AAV/plotAAVreport.R Lot-PB19045--summarize-aav-alignment

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

[1] "Lot-PB19045--summarize-aav-alignment"
Error in file(file, "rt") : invalid 'description' argument
Calls: read.table -> file
Execution halted

Thank you so much for these scripts and your help!
Warmest,
Olga

[Magdoll]

Hi Olga - please wait a few days. I'm actually fixing some code since some of the current classifications aren't correct - will update this issue when it's ready for you to test!

-Liz

[olgabot]

Hi Liz,
Thank you so much for the reply! I will wait.

While we wait, do you trust the minimap2 alignment (with the parameters from the tutorial)? It is creating a different result from what we got from the CRO and I'm wondering which data we should be more reliant on.

Thank you!
Warmest,
Olga

[Magdoll]

Hi Olga,

That is one of the things I'm playing with. I was recommending minimap2 -ax splice but now I think that is sub-optimal, minimap2 doing unnecessary "skips" (large breaks with a N cigar string) immediately with a large insertion --- in essence I think the splice mode is not working. I'm testing switching back to minimap -a instead.

[Magdoll]

Hi @olgabot - I've updated the scripts + wiki
https://github.com/Magdoll/AAV/wiki/Tutorial:-Analyzing-AAV-Data

Can you give it a go?
Thanks!
-Liz

[olgabot]

Hi @Magdoll, sure! I'll let you know if I run into any issues.

For the alignment, can you help me understand why aligning simultaneously to plasmids and hg38 is preferred? I saw that secondary alignment was turned off, which I agree with for the human alignment, but could it be useful for the plasmid alignments?

e.g. here is what I am thinking:

minimap2 --eqx -a \
   -t <threads> \
   plasmids.fasta <movie>.hifi_reads.fasta > plasmid_mapped.sam

Then get unaligned reads with:

samtools view -f 4 plasmid_mapped.sam | samtools fastq | gzip -c - > plasmid_unmapped.fastq.gz

Then map to hg38 with:

minimap2 --eqx --secondary=no -a -t 4 hg38.fa plasmid_unmapped.fastq.gz > hg38_mapped.sam

Then we can use two orthogonal methods to quantify gene impurities, and test for oncogenicity:

• Salmon/Kallisto on the plasmid unmapped fastq to quantify genes
• Featurecounts on the hg38 mapped fastq to count reads in genes

What do you think?

Thank you so much!
Warmest,
Olga

[olgabot]

Hello @Magdoll! I tested out the new code :) summarize_AAV_alignment.py worked fine for me, but I'm having issues with the report generation:

(gene-therapy)
 Wed 27 Jul - 20:35  /data/bbio-gt/2022-06-30--raav-seq-pilot/2022-07-27--test-simultaneous-alignment/03-magdoll-aav 
 olgabot@x86_64-conda-linux-gnu  Rscript ~/AAV/plotAAVreport.R $ID ../../ref/2022-jul--plasmids-plus-hg38--expr-vector-new-origin--combined-annotation.txt

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

[1] "Lot-PB19045"
`summarise()` has grouped output by 'pos0_div'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'type'. You can override using the `.groups` argument.
Warning message:
In max(type_len) : no non-missing arguments to max; returning -Inf
`summarise()` has grouped output by 'assigned_type'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'assigned_type'. You can override using the `.groups` argument.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Error: Must request at least one colour from a hue palette.
Execution halted

I can't seem to pin it down to a specific line number, but here is the error message. Thanks!

[olgabot]

Narrowed down the error to be on these lines:

  if (dim(x.read.repcap)[1] > 10) { # only plot if at least 10 reads
    grid.arrange(p1.map_starts.repcap, p1.map_ends.repcap, p1.map_len.repcap, ncol=1)
  }

Here is the error message again:

Error: Must request at least one colour from a hue palette.
Execution halted

[olgabot]

Looking at the x.reads.repcap dataframe, it seems to be that the all the assigned subtypes of repcap are NA:

333        repcap             <NA>
334        repcap             <NA>
335        repcap             <NA>
336        repcap             <NA>
337        repcap             <NA>
338        repcap             <NA>
339        repcap             <NA>
340        repcap             <NA>
341        repcap             <NA>
342        repcap             <NA>
343        repcap             <NA>

[olgabot]

For the annotation.txt, what is the "Region" for the RepCap plasmid? Is that the Rep or the Cap? And what if it overlaps the origin region, e.g. from 1000:100 in the plasmid?

NAME=myRepCap_plasmid;TYPE=repcap;REGION=500-1500;

[Magdoll]

So for repcap genome the “target” region is where the rep and cap genes are. But in all honesty I wasn’t intending to use them much. I see now where the error message is coming from and since the target region for repcap is optional, I can make sure the R script ignores some plots if subtype is all NA. Only for the vector is the “target region” required

-Liz

On Tue, Aug 2, 2022 at 10:12 AM Olga Botvinnik ***@***.***> wrote: For the annotation.txt, what is the "Region" for the RepCap plasmid? Is that the Rep or the Cap? And what if it overlaps the origin region, e.g. from 1000:100 in the plasmid? NAME=myRepCap_plasmid;TYPE=repcap;REGION=500-1500; — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAEQE356NH2GFEGK7AE3ZHLVXFJHPANCNFSM52XZYSKQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Sent from Gmail Mobile. Excuse any possible typos.

[Magdoll]

Hi Olga, please check your twitter DM as I suggested some options…

On Wed, Jul 27, 2022 at 1:44 PM Olga Botvinnik ***@***.***> wrote: Hello @Magdoll <https://github.com/Magdoll>! I tested out the new code :) summarize_AAV_alignment.py worked fine for me, but I'm having issues with the report generation: (gene-therapy) Wed 27 Jul - 20:35  /data/bbio-gt/2022-06-30--raav-seq-pilot/2022-07-27--test-simultaneous-alignment/03-magdoll-aav  ***@***.***_64-conda-linux-gnu  Rscript ~/AAV/plotAAVreport.R $ID ../../ref/2022-jul--plasmids-plus-hg38--expr-vector-new-origin--combined-annotation.txt Attaching package: ‘dplyr’ The following objects are masked from ‘package:stats’: filter, lag The following objects are masked from ‘package:base’: intersect, setdiff, setequal, union Attaching package: ‘gridExtra’ The following object is masked from ‘package:dplyr’: combine [1] "Lot-PB19045" `summarise()` has grouped output by 'pos0_div'. You can override using the `.groups` argument. `summarise()` has grouped output by 'type'. You can override using the `.groups` argument. Warning message: In max(type_len) : no non-missing arguments to max; returning -Inf `summarise()` has grouped output by 'assigned_type'. You can override using the `.groups` argument. `summarise()` has grouped output by 'assigned_type'. You can override using the `.groups` argument. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. Error: Must request at least one colour from a hue palette. Execution halted I can't seem to pin it down to a specific line number, but here is the error message. Thanks! — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAEQE357UEZ7RGRJHJG5BFTVWGNTXANCNFSM52XZYSKQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Sent from Gmail Mobile. Excuse any possible typos.