---
title: "Allele tools benchmarking"
author: "A.Mikelov"
format:
  html:
    grid: 
      body-width: 2000px
      sidebar-width: 200px
      margin-width: 200px
editor: visual

execute:
  warning: false
  message: false
  error: false
---

## Libraries, helper functions, references, colors

```{r}
#load libraries
{
  library(dplyr)
  library(tibble)
  library(readr)
  library(magrittr)
  library(tidyr)
  library(stringr)
  library(ggplot2)
  library(fuzzyjoin)
  library(plotly)
  library(seqinr)
  library(msa)
  library(conflicted)
  library(ggrepel)
  library(cowplot)
  library(tigger)
  library(bench)
  
  conflicts_prefer(dplyr::filter)
  conflicts_prefer(dplyr::first)
  conflicts_prefer(dplyr::count)
  conflicts_prefer(dplyr::select)
  conflicts_prefer(dplyr::rename)
  conflicts_prefer(base::unique)
    
  `%+%` <- function(a, b) paste0(a, b)
  
  mi_bright<-c("#42B842",
    "#845CFF",
    "#FF9429",
    "#27C2C2",
    "#E553E5",
    "#95C700",
    "#2D93FA",
    "#F05670",
    "#929BAD")
  
  mi_dark<-c(
  "#198020",
  "#5F31CC",
  "#C26A27",
  "#068A94",
  "#A324B2",
  "#659406",
  "#105BCC",
  "#AD3757",
  "#5E5E70"
  
)

mi_27<-c(
  "#99E099","#198020","#42B842","#C1ADFF","#845CFF","#5F31CC","#FFCB8F",
  "#FF9429","#C26A27","#90E0E0","#27C2C2","#068A94","#FAAAFA","#E553E5",
  "#A324B2","#CBEB67","#95C700","#659406","#99CCFF","#2D93FA","#105BCC",
  "#FFADBA","#F05670","#AD3757","#D3D7E0","#929BAD","#5E5E70"
)

##Load IUIS-named allele sequences ----
iuis_ref<-readDNAStringSet("~/ref/imgt-repseqio_100423/imgt_vj.fasta")
iuis_ref<-tibble(alleleSeq=as.character(iuis_ref),
                 alleleName=names(iuis_ref) %>% str_remove("\\|.*"))


##Load starting reference for MiXCR (assembled from C.Watson lab ref genome) ----
startingRef<-readDNAStringSet("~/.mixcr/libraries/watsonlab_lib.v.fasta")  %>% 
  as.character() %>% 
  c(readDNAStringSet("~/.mixcr/libraries/watsonlab_lib.j.fasta") %>% 
      as.character())

}
```

# Upstream processing AIRR-seq data

## Trim primers from Non-UMI data

```{bash}
ls *R1.fastq| parallel -j 6 --line-buffer 'seqtk trimfq -b 19 {} > primerTrimmed/{}'
```

## Prepare a table for downsampling raw reads

```{r}
depth_corey<-read_tsv("~/paper_alleles/oscar-coreyw_stats.tsv",
                      col_names = c("sra","donor","gz","umi_matched","totalRows","totalReads","fractionMatched"))

tibble(sra=rep(depth_corey$sra,5) %>% sort(),
       depth=rep(c(10000,25000,50000,100000,500000),length.out=360)) %>% 
  left_join(depth_corey %>% select(sra,totalReads),
            by="sra") %>%
  filter(depth<=totalReads) %>%
  select(sra,depth) %>% 
  write_tsv("~/paper_alleles/benchmarking/oscar-coreyw-downsampling_depth.tsv",col_names = F)
```

## Downsample reads with [seqtk]{.underline}

```{bash}
cat oscar-coreyw-downsampling_depth.tsv | parallel --colsep '\t' -j3 'sample_depth_k=$(echo "{2}/1000" | bc); \
	input_files=( $(find ~/paper_alleles/raw0/watson/5race -type f -name "{1}*R2*.fastq") \
	$(find ~/paper_alleles/raw0/watson/noumi/ -type f -name "{1}*R2*.fastq") ); \
	input_file="${input_files}"; \
	output_file="${input_file%%.fastq}_${sample_depth_k}K.fastq"; \
	output_file=$(echo "$output_file" | sed "s|noumi|downsampled/noumi|g"); \
		output_file=$(echo "$output_file" | sed "s|5race|downsampled/5race|g"); \
	echo $input_file;
	echo $output_file ;
	seqtk sample -s1707 "$input_file" "{2}" > "$output_file"'


cat oscar-coreyw-downsampling_depth.tsv | parallel --colsep '\t' -j3 'sample_depth_k=$(echo "{2}/1000" | bc); \
	input_files=( $(find ~/paper_alleles/raw0/watson/5race -type f -name "{1}*R1*.fastq") \
	$(find ~/paper_alleles/raw0/watson/noumi/ -type f -name "{1}*R1*.fastq") ); \
	input_file="${input_files}"; \
	output_file="${input_file%%.fastq}_${sample_depth_k}K.fastq"; \
	output_file=$(echo "$output_file" | sed "s|noumi|downsampled/noumi|g"); \
		output_file=$(echo "$output_file" | sed "s|5race|downsampled/5race|g"); \
	echo $input_file;
	echo $output_file ;
	seqtk sample -s1707 "$input_file" "{2}" > "$output_file"'	
	


	
```

## [MiXCR]{.underline}

### Upstream pipeline

**For full and downsampled data using minimal library (=1 allele per V/J gene)**

```{bash}


#Data without UMIs. "{FR1Begin(+6):FR4End}" Because the primers were on Leader and were cut
#Also adjusted quality settings in order not to take it into account because quality line
#is corrupted for these files in SRA. 
cd ~/paper_alleles/raw0/watson/noumi
ls *_R1*fastq |parallel --line-buffer -j 23 'echo {}; \
  mixcr -Xmx16g analyze generic-bcr-amplicon \
  --library watsonlab_lib \
  --species hs \
  --rna -f \
  -Massemble.cloneAssemblerParameters.badQualityThreshold=-1 \
  -Massemble.cloneAssemblerParameters.minimalQuality=1 \
  --rigid-left-alignment-boundary \
  --floating-right-alignment-boundary C    \
  --assemble-clonotypes-by "{FR1Begin(+6):FR4End}" \
  --split-clones-by C \
  {} {= s:_R1:_R2: =}\
  ~/paper_alleles/benchmarking/mixcr/noumi/{= s:_R1.fastq:: =}'

#Data without UMIs downsampled. 
cd /home/amikelov/paper_alleles/raw0/watson/downsampled/noumi

ls *_R1_*fastq |parallel --line-buffer -j 23 'echo {}; \
  mixcr -Xmx16g analyze generic-bcr-amplicon \
  --library watsonlab_lib \
  --species hs \
  --rna -f \
  -Massemble.cloneAssemblerParameters.badQualityThreshold=-1 \
  -Massemble.cloneAssemblerParameters.minimalQuality=1 \
  --rigid-left-alignment-boundary \
  --floating-right-alignment-boundary C    \
  --assemble-clonotypes-by "{FR1Begin(+6):FR4End}" \
  --split-clones-by C \
  {} {= s:_R1:_R2: =}\
  /home/amikelov/paper_alleles/benchmarking/mixcr/downsampled/noumi/{= s:_R1::;s:.fastq:: =}'

#

#Data with UMIs
cd /home/amikelov/paper_alleles/raw0/watson/5race
ls *_R1.fastq*| parallel --line-buffer -j 26 \
  'mixcr -Xmx20g analyze generic-bcr-amplicon-umi \
  -f --species hs --rna \
  --library watsonlab_lib  \
  --rigid-left-alignment-boundary \
  --floating-right-alignment-boundary C \
  --tag-pattern "^(R1:*) \ ^(UMI:N{"12"})GTACG{"1:10"}(R2:*)" \
  --assemble-clonotypes-by VDJRegion \
  --split-clones-by C \
  {} {= s:_R1:_R2: =} \
  ~/paper_alleles/benchmarking/mixcr/5race/{= s:_R1.fastq:: =}'
  
#Data without UMIs
cd /home/amikelov/paper_alleles/raw0/watson/downsampled/5race
ls *_R1_*| parallel --line-buffer -j 40 \
  'mixcr -Xmx10g analyze generic-bcr-amplicon-umi \
  -f --species hs --rna \
  --library watsonlab_lib  \
  --rigid-left-alignment-boundary \
  --floating-right-alignment-boundary C \
  --tag-pattern "^(R1:*) \ ^(UMI:N{"12"})GTACG{"1:10"}(R2:*)" \
  --assemble-clonotypes-by VDJRegion \
  --split-clones-by C \
  {} {= s:_R1_:_R2_: =} \
  ~/paper_alleles/benchmarking/mixcr/downsampled/5race/{= s:_R1::;s:.fastq:: =}'
```

### findAlleles inference with benchmarking

```{bash}

cd ~/paper_alleles/benchmarking/mixcr/noumi
mkdir findAlleles

ls *clns |parallel -j 46 --line-buffer \
  "/usr/bin/time -f '\t%e real,\t%U user,\t%S sys,\t%K avmem,\t%M maxmem' \
  -o findAlleles/{.}.time.log mixcr findAlleles \
  -t 10 \
  --no-clns-output \
  --export-library findAlleles/{.}.fasta \
  --export-alleles-mutations findAlleles/{.}.tsv \
  {} > findAlleles/{.}.report.log"
  
cd ~/paper_alleles/benchmarking/mixcr/5race
mkdir findAlleles

ls *clns |parallel -j 46 --line-buffer \
  "/usr/bin/time -f '\t%e real,\t%U user,\t%S sys,\t%K avmem,\t%M maxmem' \
  -o findAlleles/{.}.time.log mixcr findAlleles \
  -t 10 \
  --no-clns-output \
  --export-library findAlleles/{.}.fasta \
  --export-alleles-mutations findAlleles/{.}.tsv \
  {} > findAlleles/{.}.report.log"
  
cd ~/paper_alleles/benchmarking/mixcr/downsampled/noumi
mkdir findAlleles

ls *clns |parallel -j 46 --line-buffer \
  "/usr/bin/time -f '\t%e real,\t%U user,\t%S sys,\t%K avmem,\t%M maxmem' \
  -o findAlleles/{.}.time.log mixcr findAlleles -f \
  -t 10 \
  --no-clns-output \
  --export-library findAlleles/{.}.fasta \
  --export-alleles-mutations findAlleles/{.}.tsv \
  {} > findAlleles/{.}.report.log"
  
cd ~/paper_alleles/benchmarking/mixcr/downsampled/5race
mkdir findAlleles

ls SRR23446575*500K*clns |parallel -j 46 --line-buffer \
  "/usr/bin/time -f '\t%e real,\t%U user,\t%S sys,\t%K avmem,\t%M maxmem' \
  -o findAlleles/{.}.time.log mixcr -V develop findAlleles \
  -t 10 \
  --no-clns-output \
  --export-library findAlleles/{.}.fasta \
  --export-alleles-mutations findAlleles/{.}.tsv \
  {} > findAlleles/{.}.report.log"
  
```

## Immcantation upstream pipeline for full data, minimal library

### Run IgBlast

```{bash}
docker run --name igb -v /home/amikelov/:/home/amikelov/ -it immcantation/suite:4.4.0 /bin/bash
# prepare minimal reference library for igblast inside the container
mv /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta /home/amikelov/paper_alleles/
cp /home/amikelov/paper_alleles/watsonlab_IGH_VRegion_padded.fasta /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.v.fasta /usr/local/share/igblast/fasta/imgt_human_ig_v.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.j.fasta /usr/local/share/igblast/fasta/imgt_human_ig_j.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.d.fasta /usr/local/share/igblast/fasta/imgt_human_ig_d.fasta

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_v.fasta \
    -out /usr/local/share/igblast/database/imgt_human_ig_v
# D segment database

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_d.fasta \
    -out /usr/local/share/igblast/database/imgt_human_ig_d
# J segment database

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_j.fasta \
    -out /usr/localshare/igblast/database/imgt_human_ig_j

realpath  *[2s].fasta | parallel -j 88 --line-buffer  'docker exec igb changeo-igblast -r /usr/local/share/germlines/imgt/human/vdj -b /usr/local/share/igblast -s {} -n {} -p 3 -o /home/amikelov/paper_alleles/benchmarking/immcantation_fasta/igblast'
```

## Immcantation upstream downsampled data, minimal library

### Data with UMI

#### Quality filtering

```{bash}

cd /paper_alleles/raw0/watson/downsampled/5race


ls *fastq |parallel --line-buffer -j 48 \
  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ \
    immcantation/suite:4.4.0 FilterSeq.py quality \
    -s /home/amikelov/paper_alleles/raw0/watson/downsampled/5race/{} \
    -q 20 \
    --outdir /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/5race'
```

#### Extract barcodes/primers

```{bash}
ls /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/5race/*_R2_*.fastq |parallel \
  --line-buffer  -j 48 \
  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 \
             MaskPrimers.py score -s {} \
             -p /home/amikelov/paper_alleles/benchmarking/immcantation/umi.fasta \
             --start 12 --barcode --mode cut --maxerror 0.5 \
             --outname {.} \
             --log {= s::: =}.log'


/home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/5race/


realpath *_R2_*_primers-pass.fastq | parallel -j 48  --line-buffer \
  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ \
   immcantation/suite:4.4.0 PairSeq.py \
   -1 {= s:_primers-pass::;s:_R2_:_R1_: =} -2 {} \
    --2f BARCODE --coord sra'

```

#### Build consensuses

```{bash}
realpath *_R1_*pair-pass.fastq | parallel -j 48 --line-buffer 'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 BuildConsensus.py -s {} --bf BARCODE --maxerror 0.1 --outname {.} --log {.}.log'


realpath *_R2_*pair-pass.fastq | parallel -j 48 --line-buffer 'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 BuildConsensus.py -s {} --bf BARCODE --maxerror 0.1 --outname {.} --log {.}.log'


realpath *_R2_*consensus-pass.fastq | parallel -j 48 --line-buffer 'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 PairSeq.py -1 {= s:_primers-pass::;s:_R2_:_R1_: =} -2 {} --coord presto'
```

#### Assemble pairs from consensuses

```{bash}
realpath *_R2_*consensus-pass_pair-pass.fastq | parallel -j 4 --line-buffer  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ -w /home/amikelov/paper_alleles/ immcantation/suite:4.4.0 AssemblePairs.py align -1 {= s:_primers-pass::;s:_R2_:_R1_: =} -2 {} --coord presto --rc tail --1f CONSCOUNT --2f CONSCOUNT PRCONS --outname {= s:_R2::;s:_quality-pass_primers-pass_pair-pass_consensus-pass_pair-pass.fastq:: =} --log {.}.log'
```

#### Collapse UMI-duplicates

```{bash}
realpath *_assemble-pass.fastq | parallel -j 4 --line-buffer  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 ParseHeaders.py collapse -s {} -f CONSCOUNT --act min'


realpath *_assemble-pass_reheader.fastq | parallel -j 4 --line-buffer  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 CollapseSeq.py -s {} -n 20 --inner --uf PRCONS --cf CONSCOUNT --act sum --outname {.}'

```

#### Filtering to sequences with at least two representative reads

```{bash}
realpath *_collapse-unique.fastq* | parallel -j 4 --line-buffer  'docker run -v /home/amikelov/paper_alleles/:/home/amikelov/paper_alleles/ immcantation/suite:4.4.0 SplitSeq.py group -s {} -f CONSCOUNT --num 2 --outname {= s:_collapse-unique.fastq:: =}'
```

#### Run IgBlast

```{bash}
docker run --name igb -v /home/amikelov/:/home/amikelov/ -it immcantation/suite:4.4.0 /bin/bash
mv /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta /home/amikelov/paper_alleles/
cp /home/amikelov/paper_alleles/watsonlab_IGH_VRegion_padded.fasta /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.v.fasta /usr/local/share/igblast/fasta/imgt_human_ig_v.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.j.fasta /usr/local/share/igblast/fasta/imgt_human_ig_j.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.d.fasta /usr/local/share/igblast/fasta/imgt_human_ig_d.fasta

realpath *_atleast-2.fastq | parallel -j 88 --line-buffer  'docker exec igb changeo-igblast -r /usr/local/share/germlines/imgt/human/vdj -s {} -n {} -p 3 -o /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/5race/igblast'
```

### Data without UMI

#### Assemble pairs from downsampled data

```{bash}

realpath ~/paper_alleles/raw0/watson/downsampled/noumi/*_R1_*|parallel  -j 88  --line-buffer 'AssemblePairs.py align -1 {} -2 {= s:_R1_:_R2_: =} --coord sra --rc tail --outdir '~/paper_alleles/benchmarking/immcantation/downsampled/noumi/' --outname {= s:_R1::;s:.*\/:: =} --log AP.log'



```

#### Run IgBlast

```{bash}
docker run --name igb -v /home/amikelov/:/home/amikelov/ -it immcantation/suite:4.4.0 /bin/bash
mv /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta /home/amikelov/paper_alleles/
cp /home/amikelov/paper_alleles/watsonlab_IGH_VRegion_padded.fasta /usr/local/share/germlines/imgt/human/vdj/imgt_human_IGHV.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.v.fasta /usr/local/share/igblast/fasta/imgt_human_ig_v.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.j.fasta /usr/local/share/igblast/fasta/imgt_human_ig_j.fasta
cp /home/amikelov/.mixcr/libraries/watsonlab_lib.d.fasta /usr/local/share/igblast/fasta/imgt_human_ig_d.fasta

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_v.fasta \
    -out /usr/local/share/igblast/database/imgt_human_ig_v
# D segment database

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_d.fasta \
    -out /usr/local/share/igblast/database/imgt_human_ig_d
# J segment database

makeblastdb -parse_seqids -dbtype nucl -in /usr/local/share/igblast/fasta/imgt_human_ig_j.fasta \
    -out /usr/localshare/igblast/database/imgt_human_ig_j


cd /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/5race
realpath  *_collapse-unique.fastq | parallel -j 88 --line-buffer  'docker exec igb changeo-igblast -r /usr/local/share/germlines/imgt/human/vdj -b /usr/local/share/igblast -s {} -n {.} -p 3 -o /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/igblast'

cd /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/noumi
realpath  *_assemble-pass.fastq | parallel -j 88 --line-buffer  'docker exec igb changeo-igblast -r /usr/local/share/germlines/imgt/human/vdj -b /usr/local/share/igblast -s {} -n {.} -p 3 -o /home/amikelov/paper_alleles/benchmarking/immcantation/downsampled/igblast'


```

## TIgGER inference

### Load reference DBs and define function for TIgGER pipeline

```{r}

germlineDb_watsonlab<-readIgFasta("/home/amikelov/paper_alleles/watsonlab_IGH_VRegion_padded.fasta")
germlineDb_imgt<-readIgFasta("/home/amikelov/paper_alleles/imgt_human_IGHV.fasta")

tigger_pipeline<-function(filenames,out_folder,germlineDb){
  
  tryCatch( {
    
    `%+%` <- function(a, b) paste0(a, b)
    extract_id<- function(x) str_remove_all(x,".*/*.*/|\\.tsv")
    id<-extract_id(filenames)

        tigger_db <- read_tsv(filenames,show_col_types = FALSE) |> 
        mutate(v_call=str_remove(v_call,",.*"))
    
      print(filenames)
    
    novel_db <-suppressWarnings(
      suppressMessages(findNovelAlleles(tigger_db,germline_db = germlineDb,nproc=10)))
    
    allele_db<-inferGenotype(tigger_db, germline_db = germlineDb, novel=novel_db) 

    genotype<-genotypeFasta(allele_db,germline_db=germlineDb,novel=novel_db)
    writeFasta(named_sequences = genotype,out_folder %+% str_replace(id,"_airr","") %+% ".tigger.fasta")

    
    allele_db<-allele_db%>%
      mutate(countInAlleles=str_split(counts,",",simplify = F) %>% 
               lapply(as.integer) %>% 
               lapply(sum) %>% unlist(),
             total=as.integer(total),
             fractionInAlleles=countInAlleles/total)
    
    write_tsv(allele_db, out_folder %+% id %+% "_tigger.tsv")
  }, error= function(e) paste(filenames,e,sep=":" ) )
}
```

### Run TIgGER w/ benchmarking

```{r}
bench_5race_downsampled <-lapply(dir("~/paper_alleles/benchmarking/immcantation/downsampled/5race",
                  pattern="_db-pass.tsv",
                  full.names = T),
       tigger_pipeline,
       "/home/amikelov/paper_alleles/benchmarking/immcantation/tigger/",
       germlineDb = germlineDb_watsonlab)

bench_noumi_downsampled <-lapply(dir("~/paper_alleles/benchmarking/immcantation/downsampled/noumi",
                  pattern="_db-pass.tsv",
                  full.names = T),
       tigger_pipeline,
       "/home/amikelov/paper_alleles/benchmarking/immcantation/tigger/",
       germlineDb = germlineDb_watsonlab)


bench_full <-lapply(dir("~/paper_alleles/benchmarking/immcantation_fasta",
                  pattern="_db-pass.tsv",
                  full.names = T),
       tigger_pipeline,
       "/home/amikelov/paper_alleles/benchmarking/immcantation_fasta/tigger/",
       germlineDb = germlineDb_watsonlab)


sample_ids<-dir("~/paper_alleles/benchmarking/immcantation/downsampled/5race",
                  pattern="_db-pass.tsv") |> 
  str_replace("K_.*","K")

names(bench_5race_downsampled)<-sample_ids

bench_df <- lapply(sample_ids, function(id){
  b<-bench_5race_downsampled[[id]]
  print(is_tibble(b))
  if(is_tibble(b)){
    tibble(sampleId=id,
           time=b$total_time,
           memory=b$mem_alloc)
  } else {
     tibble(sampleId=id,
           time=bench::as_bench_time(NA),
           memory=bench::as_bench_bytes(NA))
  }
  
}) |> bind_rows()
```

## Time benchmarking

```{r}
read_tsv(dir("~/paper_alleles/benchmarking/mixcr/5race/findAlleles/",
             pattern = "*.time.log",
             full.names = T)
)
```

## iGenotyper - long-read based genotypes

Read iGenotyper ([https://github.com/oscarlr/IGenotyper](#0)) output - .bed haplotypes

```{r}

#list donors with repseq
donorsWithRepseq<- c(
  dir("~/paper_alleles/raw0/watson/noumi/primerTrimmed",pattern="_R1") %>%
  str_remove("_R1.fastq") %>% 
  str_remove(".*?_"),
  dir("~/paper_alleles/raw0/watson/5race",pattern="_R1") %>%
  str_remove("_R1.fastq") %>% 
  str_remove(".*?_")
)



#load iGenotyper derived haplotypes from .bed files
igenotyper_haplotype<-
  read_tsv(donorsWithRepseq %>%
             sapply(function(don){
               str_subset(dir("~/paper_alleles/benchmarking/igenotyper",
                              full.names=T),
                          don%+%".bed")}),
           col_names = c('locus',
                         "positionStart",
                         'positionEnd',
                         'geneName',
                         'haplotype',
                         'typeMatch',
                         'allele',
                         'nSeqInHaplotype'),
                               id = "fileName") %>% 
  filter(str_detect(geneName,"IGHV|IGHJ")) %>% 
  mutate(donor=str_remove_all(fileName,".*\\/|.bed"),
         geneName=factor(geneName)) %>%     
  select(geneName,donor,allele) %>% 
  unique()
```

Assign sequences to haplotypes from iGenotyper

```{r}
iuis_ref<-readDNAStringSet("~/ref/imgt_vj.fasta")
iuis_ref<-tibble(alleleSeq=as.character(iuis_ref),
                 alleleName=names(iuis_ref) %>% str_remove("\\|.*"))


startingRef<-readDNAStringSet("~/.mixcr/libraries/watsonlab_lib.v.fasta")  %>% 
  as.character() %>% 
  c(readDNAStringSet("~/.mixcr/libraries/watsonlab_lib.j.fasta") %>% 
      as.character())


igenotyper_seq<-igenotyper_haplotype %>% 
  mutate(alleleName=case_when(
    str_detect(allele,".*[0-9].*") ~ geneName %+%
      "*" %+% sprintf("%02d", as.integer(allele)), 
    str_detect(allele,"A|T|G|C") ~ "novel")
    ) %>% 
  left_join(iuis_ref,by="alleleName")  %>% 
  mutate(alleleSeq=ifelse(alleleName=="novel",
                          allele %>%  spgs::reverseComplement(case = "upper"), #iGenotyper gets novel alleles from reverse strand so need to RC
                          alleleSeq),
         inStartRef=alleleSeq %in% startingRef,            #whether sequence was in starting reference
         alleleName=ifelse(alleleName=="novel",
                           geneName %+% "*novel",
                           alleleName),
         geneType=case_when(
           str_detect(geneName,"IGHV") ~ "V",
           str_detect(geneName,"IGHD") ~ "D",
           str_detect(geneName,"IGHJ") ~ "J",
         )
  ) %>%
  select(-allele) 


#Give unique names to novel alleles within each gene----
alleleNames<-igenotyper_seq %>%
  select(-donor) %>%
  unique() %>%  
  group_by(alleleName) %>% 
  mutate(number=1:n(),
         uniqAlleleName=ifelse(str_detect(alleleName,"novel"),
                               paste(alleleName, number,sep="-"),alleleName)) %>% 
  select(alleleName,uniqAlleleName,alleleSeq)


##Assign uniq names, remove truncated genes ----
igenotyper_seq %<>% 
  #Assign unique novel allele names 
  left_join(alleleNames,by=c("alleleName","alleleSeq")) %>% 
  #remove obviously truncated V genes
  filter(!(geneType=="V" & nchar(alleleSeq)<250)) %>%  
  select(-alleleName) %>% dplyr::rename(alleleName=uniqAlleleName)

```

## Add to iGenotyper the expression information from Rep-Seq using donor-specific PacBio-derived individualized reference libraries

This is done to exclude from comparison non-expressed genes or alleles. For this we create a donor-specific reference based on long reads (iGenotyper calls)

### Write fasta files for processing with MiXCR

```{r}
#Need to filter out thrashy genes (stop-codons, no C in CDR3), for this purpose we write as fasta donor-specific reference libraries from iGenotyper haplotypes and run repseqio
##Write fasta files for processing with MiXCR ----
igenotyper_seq %>% 
  mutate(alleleName=">" %+% alleleName %+% "\n" %+% alleleSeq) %>% 
  group_by(donor,geneType) %>% 
  group_walk( ~ write_tsv(.x %>% select(alleleName),
                          "~/paper_alleles/benchmarking/igenotyper/mixcr/" %+%
                            .y$donor %+% 
                            "_iGenotyper." %+% 
                            .y$geneType %+%
                            ".fasta",
                          col_names = F))
```

### Create MiXCR formatted library with command line tool `repseqio`

```{bash}
cd ~/paper_alleles/benchmarking/igenotyper/mixcr/

ls *V*fasta|parallel -j 10 'repseqio -V alleles fromFasta -f --taxon-id 9606 --species-name hs --species-name homsap --species-name hsa --chain IGH --name-index 0 --gene-type V --gene-feature VRegion {} {.}.json'

ls *J*fasta|parallel -j 10 'repseqio -V alleles fromFasta -f --taxon-id 9606 --species-name hs --species-name homsap --species-name hsa --chain IGH --name-index 0 --gene-type J --gene-feature JRegion {} {.}.json'


```

### Merge jsons for V D J C and infer Points

```{bash}
rm -r *mifdx
cp ~/.mixcr/libraries/watson*fasta ./
cp ~/.mixcr/libraries/IGHC_repseqio.fasta ./

#Ignore "no such file" which happens for some donors - it's just because they don't have calls from PacBio for Js
ls *V*json |cut -f 1 -d "_" | parallel -j 10 --line-buffer 'echo {}; J=$(ls {}_*J.json); repseqio -V alleles  merge --replace-duplicate-genes -f  {}_*V.json /home/amikelov/.mixcr/libraries/watsonlab_lib.v.json /home/amikelov/.mixcr/libraries/watsonlab_lib.d.json $J /home/amikelov/.mixcr/libraries/watsonlab_lib.j.json /home/amikelov/.mixcr/libraries/watsonlab_lib.c.json {}_iGenotyper.json; repseqio  -V alleles  inferPoints -f -g VRegion -g JRegion  {}_iGenotyper.json {}_iGenotyper.json'
```

### Run debug and collect broken V and J genes

```{bash}

rm -r *mifdx

ls *iGenotyper.json| parallel 'repseqio -V alleles debug  -p {}' |grep -A3 IGHV > debug.txt
ls *iGenotyper.json| parallel 'repseqio -V alleles debug  -p {}' |grep -A3 IGHJ >> debug.txt

paste <(grep IGHV debug.txt|sed 's/ .*//') <(cat debug.txt | grep -E "stop|start") |sort |uniq >debug.tsv

#repseqio throws exception on  IGHV4-61*novel-2 and IGHJ1*novel-1 which are both extremely truncated iGenotyper hits. need to remove them manually
# add to debug tsv  IGHV4-61*novel-2  IGHJ1*novel-1
echo -e 'IGHV4-61*novel-2\tFR1 contains a stop codon' >> debug.tsv
echo -e 'IGHJ1*novel-1\ttruncated' >> debug.tsv
```

### Remove all non-functional/truncated genes (which didn't pass debug) and write cleaned fasta

```{r}

igenotyper_seq %<>%
  anti_join(read_tsv("~/paper_alleles/benchmarking/igenotyper/mixcr/debug.tsv",
                     col_names = c("alleleName","problem")),
            by="alleleName" )

igenotyper_seq %>% 
  mutate(alleleName=">" %+% alleleName %+% "\n" %+% alleleSeq) %>% 
  group_by(donor,geneType) %>% 
  group_walk( ~ write_tsv(.x %>% select(alleleName),
                          "~/paper_alleles/benchmarking/igenotyper/mixcr/" %+% .y$donor %+% "_iGenotyper." %+% .y$geneType %+% ".fasta",
                          col_names = F))
```

### Again create MiXCR formatted library on cleaned fasta

```{bash}
cd ~/paper_alleles/benchmarking/igenotyper/mixcr/

rm *mifdx

ls *V*fasta|parallel -j 10 'repseqio -V alleles fromFasta -f --taxon-id 9606 --species-name hs --species-name homsap --species-name hsa --chain IGH --name-index 0 --gene-type V --gene-feature VRegion {} {.}.json'

ls *J*fasta|parallel -j 10 'repseqio -V alleles fromFasta -f --taxon-id 9606 --species-name hs --species-name homsap --species-name hsa --chain IGH --name-index 0 --gene-type J --gene-feature JRegion {} {.}.json'


```

### Merge jsons for V D J C and infer Points

```{bash}
rm -r *mifdx

#Ignore "no such file" which happens for some donors - it's just because they don't have calls from PacBio for Js
ls *V*json |cut -f 1 -d "_" | parallel -j 10 --line-buffer 'echo {}; J=$(ls {}_*J.json); repseqio -V alleles  merge --replace-duplicate-genes -f  {}_*V.json /home/amikelov/.mixcr/libraries/watsonlab_lib.v.json /home/amikelov/.mixcr/libraries/watsonlab_lib.d.json $J /home/amikelov/.mixcr/libraries/watsonlab_lib.j.json /home/amikelov/.mixcr/libraries/watsonlab_lib.c.json {}_iGenotyper.json; repseqio  -V alleles  inferPoints -f -g VRegion -g JRegion  {}_iGenotyper.json {}_iGenotyper.json'
```

### Run Mixcr upstream for benchmarking data, but on donor-specific references from long-reads based genotypes (cleaned from non-functional alleles)

```{bash}

ls /home/amikelov/paper_alleles/raw0/watson/noumi/*R1.fastq | \
    cut -f 1-3 -d "_"| \
    parallel  --line-buffer -j 23 \
    'cd /home/amikelov/paper_alleles/benchmarking/igenotyper/mixcr;
     donor=$(echo {} |grep -E -o  "W-[0-9]{1,3}");
    echo {}; echo $donor;
    mixcr -Xmx22g analyze generic-bcr-amplicon -f \
    -Massemble.cloneAssemblerParameters.badQualityThreshold=1 \
    -Massemble.cloneAssemblerParameters.minimalQuality=1 \
    -Malign.parameters.vParameters.parameters.maxHits=4 \
    --species hs --rna  \
    --library ${donor}_iGenotyper \
    --rigid-left-alignment-boundary \
    --floating-right-alignment-boundary C \
    --assemble-clonotypes-by "{FR1Begin(+6):FR4End}" \
    --split-clones-by C {}*_R1.fastq {}*_R2.fastq \
    {= s:.*/::;s:_R1.fastq:: =}'


ls  /home/amikelov/paper_alleles/raw0/watson/5race/*R1*.fastq | parallel \
    --line-buffer \
    -j 27 \
    'cd /home/amikelov/paper_alleles/benchmarking/igenotyper/mixcr;
     donor=$(echo {} |grep -E -o  "S?[WC]-[0-9]{1,3}");
    echo {}; echo $donor;
    mixcr -Xmx22g analyze generic-bcr-amplicon-umi -f \
    --species hs --rna \
    --rigid-left-alignment-boundary \
    --floating-right-alignment-boundary C \
    --library ${donor}_iGenotyper \
    --tag-pattern "^(R1:*) \ ^(UMI:N{"12"})GTACG{1:10}(R2:*)"   \
    --assemble-clonotypes-by VDJRegion \
    --split-clones-by C \
    -Malign.parameters.vParameters.parameters.maxHits=4 \
    {} {= s:_R1:_R2: =} \
    {= s:.*/::;s:_R1.fastq:: =}'







```

### Read MiXCR upstream output on donor-specific references from long-reads based genotypes

```{r}

clonesets_igenotyper_files<-dir("/home/amikelov/paper_alleles/benchmarking/igenotyper/mixcr",
                                pattern="clones_IGH.tsv",
                                full.names = T)

#clonesets_igenotyper because reference is from iGenotyper

clonesets_igenotyper<-lapply(clonesets_igenotyper_files, function(fl){
  read_tsv(fl,id= "fileName") %>% 
    select(-starts_with("minQual"),-(allCAlignments:aaSeqFR4))
}) %>% bind_rows() %>% 
   mutate(donor=str_extract(fileName,"(SC|W)-[0-9]{1,3}"))

#extract mutations and donor labels
clonesets_igenotyper%<>% 
  mutate(vMut=str_split_i(allVAlignments,"\\|",i=6),
         jMut=str_split_i(allJAlignments,"\\|",i=6),
         donor=str_remove(donor,".*_"),
  )
```

### 

### Collect statistics for each V- and J-allele on \# clonotypes and \# naives

```{r}
#NOTE: I select all genes with the best alignment scores for a clone
geneStats<-clonesets_igenotyper %>% 
# select ALL best V hits and pick ones with max score.
# total # rows will be increased due to presence of very similar genes but that's ok cause these stats are only used to filter out non-expressed alleles and non-reliable allele calls 
  separate_rows(allVHitsWithScore,sep=",") %>%
  separate(allVHitsWithScore,
           into=c("vHit","score"),
           sep="\\(") %>% 
  mutate(score=as.numeric(str_remove(score,"\\)")),
         vGene=str_remove(vHit,"\\*.*")) %>% 
  group_by(fileName,cloneId) %>%
  mutate(maxScore=max(score),.after=score) %>% 
  filter(score==maxScore) %>% 
  ungroup() %>%  
#calculate nClones and nNaives for each allele  
  group_by(donor,vGene,vHit) %>% 
  summarize(nClonesInAllele=n(),
         nNaive=sum(vMut==""&jMut==""),
         cloneIds=list(cloneId)) %>% 
  ungroup() %>%
#calculate allele disbalance within each gene
  group_by(donor,vGene) %>%
  mutate(geneTotalClones= cloneIds%>% unlist() %>% unique() %>% length() ,
         alleleFreqInGene=nClonesInAllele/geneTotalClones
         ) %>%
  ungroup() %>%
  group_by(donor) %>% 
  mutate(vGeneFreq=geneTotalClones/sum(geneTotalClones)) %>% 
  ungroup() %>% 
  dplyr::rename(alleleName=vHit)

##Collect # clones and # naive clones for J genes ----
geneStats<-
  bind_rows(geneStats,
            clonesets_igenotyper %>% 
              separate_rows(allJHitsWithScore,sep=",") %>%
              separate(allJHitsWithScore,
                       into=c("jHit","score"),
                       sep="\\(") %>% 
              mutate(score=as.numeric(str_remove(score,"\\)"))) %>% 
              group_by(fileName,cloneId) %>%
              mutate(maxScore=max(score),.after=score) %>% 
              filter(score==maxScore) %>% 
              ungroup() %>%
              group_by(donor,jHit) %>% 
              summarise(nClonesInAllele=n(),
                        nNaive=sum(vMut==""&jMut==""))%>% 
              dplyr::rename(alleleName=jHit)
  )

#add it to iGenotyper sequences table
igenotyper_seq %<>%
  left_join(geneStats,by=c("donor","alleleName")) %>%
  replace_na(list(nClonesInAllele=0,nNaive=0)) 
```

## Stats TIgGER

```{r}
tiggerStats<-
  dir("~/paper_alleles/benchmarking/immcantation_fasta",
    pattern = "_db-pass.tsv",
    full.names = T) %>% 
  c(dir("~/paper_alleles/benchmarking/immcantation/downsampled/noumi",
    pattern = "_db-pass.tsv",
    full.names = T),
    dir("~/paper_alleles/benchmarking/immcantation/downsampled/5race",
    pattern = "_db-pass.tsv",
    full.names = T)) %>% 
  lapply(function(flnm){
    db<-read_tsv(flnm) %>% 
      mutate(CDR3Begin=str_locate(sequence_alignment,junction)[,1],
             CDR3VPartLen=v_germline_end-CDR3Begin) %>% 
      separate_rows(v_call,sep=",") %>% 
      mutate(geneName=str_remove(v_call,"\\*.*")) %>% group_by(geneName) %>%
      summarise(trimSeqBy=max(CDR3VPartLen),
                nClones=n()) %>%
      mutate(fileName=flnm,
             donor=str_extract(flnm,"(SC|W)-[0-9]{1,3}"),
             depth=str_extract(flnm,"10K|25K|50K|100K|500K"))
  }) %>% bind_rows() %>% 
  replace_na(list(depth="full")) %>% 
  complete(depth,nesting(donor,geneName),fill = list(nClones=0)) %>% 
  mutate(tool='tigger',
         trimSeqBy=ifelse(trimSeqBy<0,0,trimSeqBy)) 

write_tsv(tiggerStats,"~/paper_alleles/benchmarking/tiggerClonesetStatsByV.tsv")


```

## Stats MiXCR

```{r}
mixcrStats<-
  dir("~/paper_alleles/benchmarking/mixcr",
    pattern = ".clones_IGH.tsv",
    full.names = T,
    recursive = T) |> 
  lapply(function(flnm){
    db<-read_tsv(flnm,show_col_types = F) %>% 
      separate_rows(allVHitsWithScore,sep=",") %>% 
      mutate(geneName=str_remove(allVHitsWithScore,"\\*.*")) %>% group_by(geneName) %>%
      summarise(nClones=n()) %>%
      mutate(fileName=flnm,
             donor=str_extract(flnm,"(SC|W)-[0-9]{1,3}"),
             depth=str_extract(flnm,"10K|25K|50K|100K|500K"))
  }) %>% bind_rows() %>% 
  replace_na(list(depth="full")) %>% 
  complete(depth,nesting(donor,geneName),fill = list(nClones=0)) %>% 
  mutate(tool='mixcr') 

write_tsv(mixcrStats,"~/paper_alleles/benchmarking/mixcrClonesetStatsByV.tsv")
```

## Read in MIXCR [reports]{.underline} stats on clonesets

```{r}
stats<- read_tsv("~/paper_alleles/benchmarking/mixcr_stats.tsv") %>% 
  mutate(donor = str_extract(fileName,"(SC|W)-[0-9]{1,3}"))
```

## Read in MiXCR-generated statistics on allele calling

```{r}
# Read in tsv output of MiXCR findAlleles  ----
mixcr_allele_files<-c(dir("~/paper_alleles/benchmarking/mixcr/noumi/findAlleles",
                          pattern = ".tsv",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/5race/findAlleles",
                          pattern = ".tsv",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/downsampled/noumi/findAlleles",
                          pattern = ".tsv",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/downsampled/5race/findAlleles",
                          pattern = ".tsv",
                          full.names = T))

alleles_mixcr<-read_tsv(mixcr_allele_files,id="fileName") %>% 
  filter(str_detect(geneName,"IGHV|IGHJ")) %>% 
  mutate(donor = str_extract(fileName,"(SC|W)-[0-9]{1,3}"),
         depth = str_extract(fileName,"10K|25K|50K|100K|500K"),
         depth=ifelse(is.na(depth),"full",depth),
         protocol = str_extract(fileName,"noumi|5race"),
         tool="mixcr",
         trimSeqBy=str_remove_all(reliableRegion,"VRegion\\(.,-|\\)|.*\\(|\\}") %>% as.integer())
```

## Read in allele sequences from comparison tools

```{r}
fastaFromMixcr<-c(dir("~/paper_alleles/benchmarking/mixcr/noumi/findAlleles",
                          pattern = ".fasta",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/5race/findAlleles",
                          pattern = ".fasta",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/downsampled/noumi/findAlleles",
                          pattern = ".fasta",
                          full.names = T),
                     dir("~/paper_alleles/benchmarking/mixcr/downsampled/5race/findAlleles",
                          pattern = ".fasta",
                          full.names = T))
  
  
fastaFromTigger<-c(
  dir("~/paper_alleles/benchmarking/immcantation_fasta/tigger",
      pattern = "tigger.fasta",
      full.names=T),
 dir("~/paper_alleles/benchmarking/immcantation/tigger",
     pattern = "tigger.fasta",
    full.names=T))


alleleFastas<-c(fastaFromMixcr,
                fastaFromTigger)

# read in all fasta files
all_alleles<-lapply(alleleFastas,function(fasta){
  donorAlleleSeq<-readDNAStringSet(fasta)  
  donorAlleleSeq <- tibble(alleleName=names(donorAlleleSeq) %>%
                             str_remove("\\|.*"),
                           alleleSeq=as.character(donorAlleleSeq) %>% 
                             str_remove_all("\\."),
                           id = str_remove_all(fasta,".*/|.fasta|_db-pass"),
                           donor = str_extract(id, "(SC|W)-[0-9]{1,3}") %>%
                             str_remove("_.*"),
                           tool = str_extract(fasta,"tigger|mixcr"), 
                           depth = str_extract(fasta,"10K|25K|50K|100K|500K"))%>%           
    mutate(depth=ifelse(is.na(depth),"full",depth)) %>% 
    dplyr::filter(str_detect(alleleName,"IGHV|IGHJ")) %>% 
    mutate(geneName=str_remove(alleleName,"\\*.*")) %>% 
    relocate(id:geneName, .before=alleleName) %>% 
    unique()
})

protocol_donor<-tibble(
  donor=c(dir("~/paper_alleles/raw0/watson/5race",
              pattern="*_R1.fastq"),
          dir("~/paper_alleles/raw0/watson/noumi",
              pattern="*_R1.fastq")) |> 
    str_extract( "(SC|W)-[0-9]{1,3}") |>
    str_remove("_.*"),
  protocol=c(rep("5race",27),rep("noumi",46))
)

all_alleles<-all_alleles[sapply(all_alleles,nrow)!=0] %>% 
  bind_rows()  %>%  #add info on protocol and CDR3 reliable border coordinate. For TIgGER we do not count missed alleles in CDR3 parts at all
  left_join(
    alleles_mixcr %>%
      select(donor, 
             depth,
             tool,
             alleleName,
             trimSeqBy),
    by=c("donor","depth","tool","alleleName")) %>% 
  left_join(
    tiggerStats %>%
      select(donor, 
             depth,
             tool,
             geneName,
             trimSeqBy),
     by=c("donor","depth","tool","geneName")
) %>%
  rename(trimSeqBy=trimSeqBy.x) %>% 
  mutate(trimSeqBy=ifelse(is.na(trimSeqBy) & tool=="tigger",trimSeqBy.y,trimSeqBy)) %>% 
  left_join(protocol_donor,
            by="donor")
  
```

## Define genes which were not catched in long-reads

```{r}
##Filter out those genes, that are not present in iGenotyper output (not covered by longReads,not expressed, thrashy calls)----
notCatchedGenesInPullDown<-igenotyper_seq %>%
  dplyr::count(geneName,donor) %>% 
  mutate(geneName=factor(geneName,levels=alleles_mixcr %>% select(geneName) %>% pull(geneName) %>% unique())) %>% 
  na.omit() %>% 
  ungroup() %>% 
  complete(geneName,donor,fill = list(n=0) ) %>% 
  filter(n==0) %>% select(geneName,donor)

```

## Compare alleles by sequence

Requiring exact sequence match of substring covered in AIRR-seq.

```{r}
tools_depths<-all_alleles %>% 
  select(donor,tool,depth) %>% unique()

compareJs<-F

alleleComparison<-
  apply(tools_depths, 1, function(tldpt){
    regex_full_join(
      igenotyper_seq %>%
                      filter(donor == tldpt["donor"]) %>% 
                      filter(!str_detect(alleleName,"IGHJ")|compareJs),
      
      all_alleles %>%
                      filter(donor == tldpt["donor"],
                             tool == tldpt["tool"],
                             depth == tldpt["depth"]) %>%
                      filter(!str_detect(alleleName,"IGHJ")|compareJs) %>% 
                      mutate(alleleSeq=ifelse(is.na(trimSeqBy),
                                              alleleSeq,
                                              str_sub(alleleSeq,1,-trimSeqBy-1))),
      by =c ("alleleSeq","donor")) %>%  
      replace_na(list(donor.y = tldpt["donor"],
                      tool = tldpt["tool"],
                      depth = tldpt["depth"]))
  }) %>% bind_rows() %>% 
  mutate(geneName=ifelse(is.na(geneName.y),as.character(geneName.x),geneName.y),.before=alleleSeq.x) %>% 
  ungroup() %>%
  dplyr::rename(
    donor=donor.y,
    alleleName.iGenotyper=alleleName.x,
    alleleSeq.iGenotyper=alleleSeq.x,
    alleleName.tool=alleleName.y,
    alleleSeq.tool=alleleSeq.y) %>% 
  select(-geneName.x,-geneName.y,-donor.x,-protocol)

alleleComparison_bckp<-alleleComparison


#Add information from clonsetsfindAlleles for each gene/donor
alleleComparison %<>% 
  left_join(alleles_mixcr %>%
              select(geneName,
                     depth,
                     donor,
                     enoughInfo) %>%
              unique(),
            by= c("donor",
                  "depth",
                  "geneName")) %>%
  left_join(
    bind_rows(
      mixcrStats,
      tiggerStats
    ) %>%
      select(-fileName,-trimSeqBy) %>% 
      rename(nClonesTotalInGene=nClones),
    by=c('donor','depth','tool','geneName')
  ) %>% 
  mutate(inStartRef=ifelse(is.na(inStartRef),alleleName.tool %in% names(startingRef),inStartRef)) %>%
  relocate(tool,donor,geneType,geneName,
           alleleName.iGenotyper,alleleName.tool,
           id,depth,inStartRef,nClonesInAllele,nNaive,
           alleleSeq.iGenotyper,alleleSeq.tool,
           nClonesTotalInGene) %>%
  #add info on protocol
  left_join(protocol_donor,
            by="donor")


alleleComparison %<>%
  left_join(notCatchedGenesInPullDown %>% mutate(geneNotCatchedInPacBio=TRUE),
            by=c("geneName","donor")) %>% 
  replace_na(list(geneNotCatchedInPacBio=FALSE)) %>% 
  mutate(matchStatus=
           case_when(
             nClonesInAllele<=10 ~ "lowExpression",
             nNaive<=2 ~ "noSuchNaiveExpress",
             is.na(alleleName.iGenotyper) & geneNotCatchedInPacBio ~ "geneNotCatchedInPacBio",
             tool=="mixcr" & enoughInfo == F ~ "notEnoughInfo",
             tool=="tigger" & nClonesTotalInGene<50 & inStartRef=="FALSE" ~ "notEnoughInfo",
            !is.na(alleleName.tool) & !is.na(alleleName.iGenotyper) ~ "match",
             is.na(alleleName.iGenotyper) ~ "falsePositive",
             is.na(alleleName.tool) ~ "falseNegative"
           ))  %>% 
  #join MiXCR stats and filter out shallow sample (<500K reads)
  left_join(stats %>% select(donor,totalReads),by="donor") %>%
  filter(totalReads>500000)
```

### Genes not catched in PacBio but alleles matched to other genes

```{r}
fpSeq<-alleleComparison %>% 
  filter(tool=="mixcr",
         matchStatus=="falsePositive",
         depth=='full') %>% 
  select(alleleSeq.tool,alleleName.tool,donor)

notCatchedPBAllelesToDifGenes<-alleleComparison %>%
  filter(depth=='full') %>% 
  mutate(geneName.iGenotyper=str_remove(alleleName.iGenotyper,"\\*.*"),
         geneName.tool=str_remove(alleleName.tool,"\\*.*")) %>%
  regex_inner_join(fpSeq,by = join_by(alleleSeq.iGenotyper==alleleSeq.tool)) %>%
  filter(geneName.iGenotyper!=geneName.tool) %>% 
  select(geneName.iGenotyper,geneName.tool,donor.y) %>% 
  unique() %>% 
  rename(donor=donor.y,geneName=geneName.iGenotyper) %>%
  inner_join(notCatchedGenesInPullDown,by=c("donor","geneName")) %>% 
  select(-geneName) %>% 
  unique() %>% 
  mutate(alleleMatchedToOtherGene=TRUE)

alleleComparison %<>%
  left_join(notCatchedPBAllelesToDifGenes %>%
              mutate(alleleMatchedToOtherGene=TRUE), 
            by= join_by(donor == donor , geneName==geneName.tool)) %>% 
  replace_na(list(alleleMatchedToOtherGene=F)) %>% 
  mutate(matchStatus=ifelse(alleleMatchedToOtherGene,
                            "geneNotCatchedInPacBio",
                            matchStatus))

```

## Fig 1 a/b Comparison NOT divided by presence/absence in starting reference library

```{r}
comparisonStats<-
    alleleComparison %>%
   filter(depth=="full") %>%
  ungroup() %>% 
  group_by(donor,depth,tool,matchStatus) %>% 
  # n_disctinct because some of the tools' findings match to multiple genes, while may be the same in sequence
  summarise(N=n_distinct(alleleSeq.tool,alleleName.iGenotyper)) %>% 
  ungroup() %>% 
  complete(matchStatus,depth,nesting(donor,tool),fill=list(N=0)) %>% 
  group_by(donor,tool) %>% 
  mutate(nAllelesTotal=sum(N[matchStatus == "match" | matchStatus == "falseNegative"])) %>% 
  mutate(fraction=N/nAllelesTotal) %>% 
  left_join(stats %>% select(donor,totalReads:readsUsedInClonotypes)) |> 
  replace_na(list(fraction=0))


# Fraction matched

give.n <- function(y){
  return(data.frame(
    y = -0.07,
    label = paste0('N=', length(y)))) 
} 

give.median <- function(y){
  return(data.frame(
    y = -0.14,
    label =  paste0('med=', round(median(y),2)))) 
}

fig1a<-comparisonStats  %>% 
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("match")) %>%   
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
  ggplot(aes(y=fraction,x=tool,fill=tool))+
  geom_boxplot(outlier.shape = NA)+
  facet_grid( ~ matchStatus,
              labeller = labeller(inStartRef=
                                    c(`FALSE`="Absent in starting reference library",
                                      `TRUE`="Present in starting reference library"),
                                  matchStatus=c(match="Alleles supported\nby long reads",
                                                falseNegative="Alleles missed",
                                                falsePositive="Alleles not supported\nby long reads"))
  )+
  geom_point(position = position_dodge(width=0.75))+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="Fraction of detected alleles",
       x="Genotyping tool",
       fill="Genotyping\ntool")+
  stat_summary(fun.data = give.n, geom = "text",size=3)+
  stat_summary(fun.data = give.median, geom = "text",size=3)+
  scale_x_discrete(labels=c("MiXCR","TIgGER"))+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"),
        axis.title.x =element_blank()
        )



# N False-positives


give.n <- function(y){
  return(data.frame(
    y = -1,
    label = paste0('N=', length(y)))) 
} 

give.median <- function(y){
  return(data.frame(
    y = -2.2,
    label =  paste0('med=', round(median(y),2)))) 
}

fig1b<-comparisonStats  %>% 
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("falsePositive")) %>%  
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
  ggplot(aes(y=N,x=tool,fill=tool))+
  geom_boxplot(outlier.shape = NA)+
  facet_grid( ~ matchStatus,
              labeller = labeller(inStartRef=
                                    c(`FALSE`="Absent in starting reference library",
                                      `TRUE`="Present in starting reference library"),
                                  matchStatus=c(match="Alleles supported\nby long reads",
                                                falseNegative="Alleles missed",
                                                falsePositive="Alleles not supported\nby long reads"))
  )+
  geom_point(position = position_dodge(width=0.75))+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="# potential false positive calls",
       x="Genotyping tool",
       fill="Genotyping\ntool")+
  stat_summary(fun.data = give.n, geom = "text",size=3)+
  stat_summary(fun.data = give.median, geom = "text",size=3)+
  scale_x_discrete(labels=c("MiXCR","TIgGER"))+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"),
        axis.title.x=element_blank())



```

## Fig 1c. Comparison divided by presence/absence in starting reference library

```{r}
comparisonStats<-
  alleleComparison %>%
    filter(depth=="full") %>%
  ungroup() %>% 
  group_by(donor,inStartRef,tool,matchStatus) %>% 
  # n_disctinct because some of the tools' findings match to multiple genes, while may be the same in sequence
  summarise(N=n_distinct(alleleSeq.tool,alleleName.iGenotyper)) %>% 
  ungroup() %>% 
  complete(inStartRef,matchStatus,nesting(donor,tool),fill=list(N=0)) %>% 
  group_by(donor,tool) %>% 
  mutate(nAllelesTotal=sum(N[matchStatus=="match" | matchStatus== "falseNegative"])) %>% 
  mutate(fraction=N/nAllelesTotal)


#Comparison divided by presence/absence in starting reference library ----

give.n <- function(y){
  return(data.frame(
    y = -0.03,
    label = paste0('N=', length(y)))) 
}

give.median <- function(y){
  return(data.frame(
    y = -0.06,
    label =  paste0('median=', round(median(y),2)))) 
}

fig1c<-  comparisonStats  %>% 
    left_join(stats %>% select(donor,totalReads:readsUsedInClonotypes)) %>% 
    filter(
      totalClonotypes>3000,
      matchStatus %in% c("falseNegative")) %>% 
    mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
    ggplot(aes(y=fraction,x=tool,fill=tool))+
    geom_boxplot(outlier.shape = NA)+
    facet_grid( ~ inStartRef,
               labeller = labeller(inStartRef=c(`FALSE`="Absent in starting\nreference library",
                                                `TRUE`="Present in starting\nreference library"),
                                   matchStatus=c(match="Alleles supported by long reads",
                                                 falseNegative="Alleles missed",
                                                 falsePositive="Alleles not supported by long reads"))
               )+
    geom_point(aes(group=inStartRef),position = position_dodge(width=0.75))+
    theme_bw()+
    scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
    labs(y="Fraction of missed alleles",
         x="Genotyping tool",
         fill="Genotyping\ntool")+
    stat_summary(fun.data = give.n, geom = "text",size=3)+
    stat_summary(fun.data = give.median, geom = "text",size=3) +
    scale_x_discrete(labels=c("MiXCR","TIgGER"))+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"),
        axis.title.x=element_blank())
```

# FigS4. Comparison divided by TYPE (5'RACE/ Vmix) of data

```{r}
comparisonStats<-
  alleleComparison %>%
    filter(depth=="full") %>%
  ungroup() %>% 
  group_by(donor,inStartRef,tool,matchStatus) %>% 
  # n_disctinct because some of the tools' findings match to multiple genes, while may be the same in sequence
  summarise(N=n_distinct(alleleSeq.tool,alleleName.iGenotyper)) %>% 
  ungroup() %>% 
  complete(inStartRef,matchStatus,nesting(donor,tool),fill=list(N=0)) %>% 
  group_by(donor,tool) %>% 
  mutate(nAllelesTotal=sum(N[matchStatus=="match" | matchStatus== "falseNegative"])) %>% 
  mutate(fraction=N/nAllelesTotal)


#Comparison divided by presence/absence in starting reference library ----

give.n <- function(y){
  return(data.frame(
    y = -0.03,
    label = paste0('N=', length(y)))) 
}

give.median <- function(y){
  return(data.frame(
    y = -0.06,
    label =  paste0('median=', round(median(y),2)))) 
}

figS4 <-  comparisonStats  %>% 
    left_join(stats %>% select(donor,totalReads:readsUsedInClonotypes)) %>% 
    filter(
      totalClonotypes>3000,
      matchStatus %in% c("falseNegative")) %>% 
    mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
    ggplot(aes(y=fraction,x=tool,fill=tool))+
    geom_boxplot(outlier.shape = NA)+
    facet_grid( ~ inStartRef,
               labeller = labeller(inStartRef=c(`FALSE`="Absent in starting\nreference library",
                                                `TRUE`="Present in starting\nreference library"),
                                   matchStatus=c(match="Alleles supported by long reads",
                                                 falseNegative="Alleles missed",
                                                 falsePositive="Alleles not supported by long reads"))
               )+
    geom_point(aes(group=inStartRef),position = position_dodge(width=0.75))+
    theme_bw()+
    scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
    labs(y="Fraction of missed alleles",
         x="Genotyping tool",
         fill="Genotyping\ntool")+
    stat_summary(fun.data = give.n, geom = "text",size=3)+
    stat_summary(fun.data = give.median, geom = "text",size=3) +
    scale_x_discrete(labels=c("MiXCR","TIgGER"))+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"),
        axis.title.x=element_blank())
```

## Fig1d/e. Downsampling effects

```{r}
comparisonStats<-
  alleleComparison %>%
   ungroup() %>% 
  group_by(depth,donor,tool,matchStatus) %>% 
  summarise(N=n_distinct(alleleSeq.tool,alleleName.iGenotyper)) %>% 
  ungroup() %>% 
  mutate(
    donor= factor(donor,levels= unique(donor)),
    depth=factor(depth,levels=c("full", "500K","100K","50K","25K","10K")),
    tool=factor(tool,levels=c('mixcr','tigger'))
  ) %>% 
  complete(matchStatus,depth,donor,tool,fill=list(N=0)) %>% 
  group_by(donor,tool,depth) %>% 
  mutate(nAllelesTotal=sum(N[matchStatus=="match" | matchStatus== "falseNegative"])) %>% 
  mutate(fraction=N/nAllelesTotal) %>% 
  left_join(stats %>% select(donor,totalReads:readsUsedInClonotypes)) %>%
    left_join(protocol_donor) 

give.n <- function(y){
  return(data.frame(
    y = -0.05,
    label = paste0('N=', length(y)))) 
}
give.median <- function(y){
  return(data.frame(
    y = -0.13,
    label =  paste0('med=', round(median(y),2)))) 
}

fig1d<-comparisonStats  %>% 
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("match")) %>% 
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive")),
         depth=factor(depth,levels=c("full", "500K","100K","50K", "25K","10K"))) %>% 
  ggplot(aes(y=fraction,x=depth,fill=tool))+
  geom_boxplot(outlier.shape = NA,position = position_dodge2(preserve = "single"))+
  facet_wrap(  ~ matchStatus,
               labeller = labeller(inStartRef=c(`FALSE`="Absent in starting reference library",
                                                `TRUE`="Present in starting reference library"),
                                   matchStatus=c(match="Alleles supported by long reads",
                                                 falseNegative="Alleles missed",
                                                 falsePositive="Alleles not supported by long reads")),
               nrow=2
  )+
  geom_point(position = position_dodge(0.75, preserve = 'total'))+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="Fraction of detected alleles",
       x="",
       fill="Genotyping tool") +
  theme(legend.position = "top",
        strip.background = element_rect(fill="white"))+
  stat_summary(fun.data = give.n, geom = "text",size=3,position =position_dodge(0.75, preserve = 'total') )+
  stat_summary(fun.data = give.median, geom = "text",size=3,position =position_dodge(0.75, preserve = 'total'))


give.n <- function(y){
  return(data.frame(
    y = -1,
    label = paste0('N=', length(y)))) 
}
give.median <- function(y){
  return(data.frame(
    y = -2,
    label =  paste0('med=', round(median(y),2)))) 
}


fig1e<-comparisonStats  %>% 
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("falsePositive")) %>% 
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive")),
         depth=factor(depth,levels=c("full", "500K","100K","50K", "25K","10K"))) %>% 
  ggplot(aes(y=N,x=depth,fill=tool))+
  geom_boxplot(outlier.shape = NA,position = position_dodge2(preserve = "single"))+
  facet_wrap(  ~ matchStatus,
               labeller = labeller(inStartRef=c(`FALSE`="Absent in starting reference library",
                                                `TRUE`="Present in starting reference library"),
                                   matchStatus=c(match="Alleles supported by long reads",
                                                 falseNegative="Alleles missed",
                                                 falsePositive="Alleles not supported by long reads")),
               nrow=2
  )+
  ggbeeswarm::geom_quasirandom(dodge.width=.8)+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="# potential false positive calls",
       x="Downsampling depth",
       fill="Genotyping tool") +
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"))+
  stat_summary(fun.data = give.n, geom = "text",size=3,position =position_dodge(0.75, preserve = 'total') )+
  stat_summary(fun.data = give.median, geom = "text",size=3,position =position_dodge(0.75, preserve = 'total'))


```

## Fig1F. Detection vs allelic imbalance / usage

```{r}
alleleComparisonFull17<-
alleleComparison %>%  
  filter(matchStatus!="geneNotCatchedInPacBio") %>% 
  mutate(
         tool=case_when(tool=="mixcr" ~ "MiXCR",
                        tool=="tigger" ~ "TIgGER"),
         detected= matchStatus=="match",
         includedInAnalysis = matchStatus %in% c("match","falsePositive","falseNegative"),
         depth=factor(depth,levels=c("full","500K","100K","50K","10K")),
         ) %>% 
  filter(donor=="SC-17",!is.na(vGeneFreq)) 

fig1f<-alleleComparisonFull17 %>%  
  ggplot()+
  geom_point(aes(x=vGeneFreq,
                 y=alleleFreqInGene,
                 color=detected,
                 alpha=includedInAnalysis,
                 shape=includedInAnalysis),
             size=5)+
  scale_alpha_manual(values=c(0.3,0.75),guide="none")+
  theme_bw()+
  scale_color_manual(values = mi_bright[c(3,4)],labels=c("Not detected","Detected"))+
  scale_shape_manual(values=c(13,16),labels=c("Not enough clones, excluded","Enough clones, included"))+
  facet_wrap(tool~depth,nrow=2)+
  labs(x= "V gene frequency within the repertoire (V usage)",
       y= "Allele frequency within one V gene (Allelic Imbalance)",
       shape = "Allele expression",
       color ="Allele detection")+
  geom_text(data=alleleComparisonFull17 %>% 
              group_by(depth,tool) %>% 
              summarise(nAllelesTotal=sum(includedInAnalysis),
                        nDetected=sum(matchStatus=="match")),
              aes(label="Detected: " %+% nDetected   ),
            x = .06, y = 0.8, hjust=1, vjust=0 , size=3
  ) +
  theme(legend.position = "bottom",
        legend.box = "vertical")+ 
  scale_x_continuous(breaks = seq(0,0.06, by = 0.02), 
      labels = c(".0",".02",".04",".06"))

 
```

## Figure 1

```{r}
fig1<-plot_grid( plot_grid(fig1a,fig1b, fig1c,nrow=1,labels = c("A","B","C"),
                           rel_widths = c(1,1,2)),
                    fig1d,
                    fig1e,
                    fig1f,
                    nrow=4,
                    labels=c("","D","E","F"),
                    rel_heights = c(1,1,1,2))
ggsave("~/paper_alleles/analysis/plots/fig1.pdf",fig1,units = "in",width = 10,height = 16,bg="white")
ggsave("~/paper_alleles/analysis/plots/fig1.png",fig1,units = "in",width =10,height =16,bg="white")
```

## 

## COMPARE ON FULL REF

Same upstream, but removing \--library watsonlab_lib for MiXCR, and with standard IgBlast reference library.

### Fig S1. Read in allele sequences from comparison tools

```{r}

alleles_mixcr_full<-read_tsv(
  dir("/home/amikelov/paper_alleles/benchmarking/full_mixcr",
      pattern="tsv",
      full.names=T) %>% str_subset("clones",negate=T),
  id="fileName") %>% 
  filter(str_detect(geneName,"IGHV|IGHJ")) %>% 
  mutate( donor = str_extract(fileName,"(SC|W)-[0-9]{1,3}")) %>% 
  #should give 'Missing pieces' warning
  separate(donor,sep="_", into=c("donor","depth")) %>% 
  mutate(depth=ifelse(is.na(depth),"full",depth),
         protocol = str_extract(fileName,"noumi|5race"),
         tool='mixcr',
         trimSeqBy=str_remove_all(reliableRegion,"VRegion\\(.,-|\\)|.*\\(|\\}|-") %>% as.integer())

fastaFromMixcr_full<-
  dir("~/paper_alleles/benchmarking/full_mixcr",
      pattern="mixcr.fasta",
      full.names = T)
  
  
fastaFromTigger_imgt<-dir("~/paper_alleles/benchmarking/imgt_immcantation/tigger",
      pattern = "tigger.fasta",
      full.names=T)



alleleFastas_full<-c(fastaFromMixcr_full,
                fastaFromTigger_imgt)



# read in all fasta files
all_alleles_full<-lapply(alleleFastas_full,function(fasta){
  donorAlleleSeq<-readDNAStringSet(fasta)  
  donorAlleleSeq <- tibble(alleleName=names(donorAlleleSeq) %>%
                             str_remove("\\|.*"),
                           alleleSeq=as.character(donorAlleleSeq) %>% 
                             str_remove_all("\\."),
                           id = str_remove_all(fasta,".*/|.fasta|_db-pass"),
                           donor=str_extract(id, "(SC|W)-[0-9]{1,3}") %>%
                             str_remove("_.*"),
                           tool=str_remove(id,".*\\."),
                           depth=str_extract(fasta,"10K|50K|100K|500K"))%>%           
    mutate(depth=ifelse(is.na(depth),"full",depth)) %>% 
    dplyr::filter(str_detect(alleleName,"IGHV|IGHJ")) %>% 
    mutate(geneName=str_remove(alleleName,"\\*.*")) %>% 
    relocate(id:geneName, .before=alleleName) %>% 
    unique()
})


all_alleles_full<-all_alleles_full[sapply(all_alleles_full,nrow)!=0] %>% 
  bind_rows()  %>%  #add info on protocol and CDR3 reliable border
  left_join(
    alleles_mixcr_full %>%
      select(donor, 
             depth,
             tool,
             alleleName,
             trimSeqBy),
    by=c("donor","depth","tool","alleleName")) %>% 
  left_join(
    tiggerStats %>%
      select(donor, 
             depth,
             tool,
             geneName,
             trimSeqBy),
     by=c("donor","depth","tool","geneName")
) %>%
  rename(trimSeqBy=trimSeqBy.x) %>% 
  mutate(trimSeqBy=ifelse(is.na(trimSeqBy) & tool=="tigger",trimSeqBy.y,trimSeqBy)) %>% 
  left_join(read_tsv("~/paper_alleles/benchmarking/raw/repseq_meta.tsv") %>% 
            select(donor,protocol),
            by="donor")
  
```

### Compare alleles by sequence

```{r}
tools_depths_full<-all_alleles_full %>% 
  select(donor,tool,depth) %>% unique()

compareJs<-F

alleleComparison_full<-
  apply(tools_depths_full, 1, function(tldpt){
    regex_full_join(
      
      igenotyper_seq %>%
                      filter(donor == tldpt["donor"]) %>% 
                      filter(!str_detect(alleleName,"IGHJ")|compareJs),
      
      all_alleles_full %>%
                      filter(donor == tldpt["donor"],
                             tool == tldpt["tool"],
                             depth == tldpt["depth"]) %>%
                      filter(!str_detect(alleleName,"IGHJ")|compareJs) %>% 
                      mutate(alleleSeq=ifelse(is.na(trimSeqBy),
                                              alleleSeq,
                                              str_sub(alleleSeq,1,-trimSeqBy-1)),
                             alleleSeq=ifelse(protocol=="noumi",
                                              str_sub(alleleSeq,20),
                                              alleleSeq)),
      by =c ("alleleSeq","donor")) %>%  
      replace_na(list(donor.y = tldpt["donor"],
                      tool = tldpt["tool"],
                      depth = tldpt["depth"]))
  }) %>% bind_rows() %>% 
  mutate(geneName=ifelse(is.na(geneName.y),geneName.x,geneName.y),.before=alleleSeq.x) %>% 
  ungroup() %>%
  dplyr::rename(
    donor=donor.y,
    alleleName.iGenotyper=alleleName.x,
    alleleSeq.iGenotyper=alleleSeq.x,
    alleleName.tool=alleleName.y,
    alleleSeq.tool=alleleSeq.y) %>% 
  select(-geneName.x,-geneName.y,-donor.x,-protocol)



#Add information from findAlleles for each gene/donor
alleleComparison_full %<>% 
   left_join(alleles_mixcr_full %>%
             # filter(depth=="full") %>%
              select(geneName,
                     depth,
                     donor,
                     status,
                     enoughInfo) %>%
              unique(),
            by= c("donor",
                  "depth",
                  "geneName")) %>%

  left_join(
    bind_rows(
      mixcrStats,
      tiggerStats
    ) %>%
      select(-fileName,-trimSeqBy) %>% 
      rename(nClonesTotalInGene=nClones),
    by=c('donor','depth','tool','geneName')
  ) %>% 
  mutate(inStartRef=ifelse(is.na(inStartRef),alleleName.tool %in% names(startingRef),inStartRef)) %>%
  relocate(tool,donor,geneType,geneName,
           alleleName.iGenotyper,alleleName.tool,
           id,depth,inStartRef,nClonesInAllele,nNaive,
           alleleSeq.iGenotyper,alleleSeq.tool,
           enoughInfo) %>% 
  #add info on protocol
  left_join(read_tsv("~/paper_alleles/benchmarking/raw/repseq_meta.tsv") %>% 
              select(donor,protocol),
            by="donor")




alleleComparison_full %<>%
  left_join(notCatchedGenesInPullDown %>% mutate(geneNotCatchedInPacBio=TRUE),
            by=c("geneName","donor")) %>% 
  replace_na(list(geneNotCatchedInPacBio=FALSE)) %>% 
  mutate(matchStatus=
           case_when(
             nClonesInAllele<=10 ~ "lowExpression",
             nNaive<=2 ~ "noSuchNaiveExpress",
             is.na(alleleName.iGenotyper) & geneNotCatchedInPacBio ~ "geneNotCatchedInPacBio",
              tool=="mixcr" & enoughInfo == F & (status == "COULD_NOT_BE_ALIGNED_ON_KNOWN_VARIANT" | status == "NO_CLONES_TO_SEARCH") ~ "notEnoughInfo",
             tool=="tigger" & nClonesTotalInGene<50 & inStartRef=="FALSE" ~ "notEnoughInfo",
             !is.na(alleleName.tool) & !is.na(alleleName.iGenotyper) ~ "match",
             is.na(alleleName.iGenotyper) ~ "falsePositive",
             is.na(alleleName.tool) ~ "falseNegative"
           )) %>% 
  #join MiXCR stats and filter out shallow sample (<500K reads)
  left_join(stats %>% select(donor,totalReads),by="donor") %>%
  filter(totalReads>500000)


```

### Genes not catched in PacBio but alleles matched to other genes

```{r}
fpSeq<-alleleComparison_full %>% 
  filter(tool=="mixcr",matchStatus=="falsePositive") %>% 
  select(alleleSeq.tool,alleleName.tool,donor)

notCatchedPBAllelesToDifGenes_full<-alleleComparison_full %>%
  mutate(geneName.iGenotyper=str_remove(alleleName.iGenotyper,"\\*.*"),
         geneName.tool=str_remove(alleleName.tool,"\\*.*")) %>%
  regex_inner_join(fpSeq,by = join_by(alleleSeq.iGenotyper==alleleSeq.tool)) %>%
  filter(geneName.iGenotyper!=geneName.tool) %>% 
  select(geneName.iGenotyper,geneName.tool,donor.y) %>% 
  unique() %>% 
  rename(donor=donor.y,geneName=geneName.iGenotyper) %>%
  inner_join(notCatchedGenesInPullDown,by=c("donor","geneName")) %>% 
  select(-geneName) %>% 
  unique() %>% 
  mutate(alleleMatchedToOtherGene=TRUE)

alleleComparison_full %<>%
  left_join(notCatchedPBAllelesToDifGenes_full %>%
              mutate(alleleMatchedToOtherGene=TRUE), 
            by= join_by(donor == donor , geneName==geneName.tool)) %>% 
  replace_na(list(alleleMatchedToOtherGene=F)) %>% 
  mutate(matchStatus=ifelse(alleleMatchedToOtherGene,
                            "geneNotCatchedInPacBio",
                            matchStatus))

```

## FigS1a/b. Comparison NOT divided by presence/absence in starting reference library

```{r}
comparisonStats_full<-
    alleleComparison_full %>%
   filter(depth=="full") %>%
  ungroup() %>% 
  group_by(donor,tool,matchStatus) %>% 
  # n_disctinct because some of the tools' findings match to multiple genes, while may be the same in sequence
  summarise(N=n_distinct(alleleSeq.tool,alleleName.iGenotyper),
            geneNames=paste(unique(geneName),collapse=",")) %>% 
  ungroup() %>% 
  complete(matchStatus,nesting(donor,tool),fill=list(N=0)) %>% 
  group_by(donor,tool) %>% 
  mutate(nAllelesTotal=sum(N[matchStatus=="match" | matchStatus== "falseNegative"])) %>% 
  mutate(fraction=N/nAllelesTotal) %>% 
  left_join(stats %>% select(donor,totalReads:readsUsedInClonotypes)) 


# Fraction matched

give.n <- function(y){
  return(data.frame(
    y = -0.07,
    label = paste0('N=', length(y)))) 
} 

give.median <- function(y){
  return(data.frame(
    y = -0.14,
    label =  paste0('med=', round(median(y),4)))) 
}

figS1a<-comparisonStats_full  %>%
  mutate(
    grp=tool,
    tool = ifelse(tool=="mixcr","MiXCR Full Lib",tool),
    tool = ifelse(tool=="tigger","TIgGER Full Lib",tool)) %>% 
  bind_rows(comparisonStats %>% 
              ungroup() %>% 
              filter(depth=="full") %>% 
              mutate(grp=tool,
                     tool = ifelse(grp=="mixcr","MiXCR",tool), 
                     tool = ifelse(grp=="tigger","TIgGER",tool)) 
              ) %>% 
  left_join(read_tsv("~/paper_alleles/benchmarking/raw/repseq_meta.tsv") %>%
  select(donor,protocol),
  by="donor") %>%
  filter(
    totalClonotypes>3000,
    matchStatus == "match") %>%   
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
  ggplot(aes(y=fraction,x=tool,fill=grp))+
  geom_boxplot(outlier.shape = NA)+
  facet_grid( ~ grp,
              labeller = labeller(grp =c (mixcr="MiXCR",tigger="TIgGER")),
              scales = "free_x"
  )+
  geom_point(position = position_dodge(width=0.75))+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="Fraction of detected alleles",
       x="Genotyping tool",
       fill="Genotyping\ntool")+
  stat_summary(fun.data = give.n, geom = "text",size=3)+
  stat_summary(fun.data = give.median, geom = "text",size=3)+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"))+
  geom_line(aes(group=donor),alpha=0.5,color="darkgrey")



# N False-positives


give.n <- function(y){
  return(data.frame(
    y = -1,
    label = paste0('N=', length(y)))) 
} 

give.median <- function(y){
  return(data.frame(
    y = -2.2,
    label =  paste0('med=', round(median(y),2)))) 
}


figS1b<-comparisonStats_full  %>%
  mutate(
    grp=tool,
    tool = ifelse(tool=="mixcr","MiXCR Full Lib",tool),
    tool = ifelse(tool=="tigger","TIgGER Full Lib",tool)) %>% 
  bind_rows(comparisonStats %>% 
              ungroup() %>% 
              filter(depth=="full") %>% 
              mutate(grp=tool,
                     tool = ifelse(grp=="mixcr","MiXCR",tool),
                     tool = ifelse(grp=="tigger","TIgGER",tool))
              ) %>% 
  left_join(read_tsv("~/paper_alleles/benchmarking/raw/repseq_meta.tsv") %>%
  select(donor,protocol),
  by="donor") %>%
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("falsePositive")) %>%  
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
  ggplot(aes(y=N,x=tool,fill=grp))+
  geom_boxplot(outlier.shape = NA)+
  facet_grid( ~ grp,
              labeller = labeller(grp =c (mixcr="MiXCR",tigger="TIgGER")),
              scales = "free_x"
  )+
  geom_point(position = position_dodge(width=0.75))+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  labs(y="# potential false positive allele calls",
       x="Genotyping tool",
       fill="Genotyping\ntool")+
  stat_summary(fun.data = give.n, geom = "text",size=3)+
  stat_summary(fun.data = give.median, geom = "text",size=3)+
  theme(legend.position = "none",
        strip.background = element_rect(fill="white"))+
  geom_line(aes(group=donor),alpha=0.5,color="darkgrey")

figS1<-cowplot::plot_grid(figS1a,figS1b,nrow=2,labels=c("A","B"))
ggsave("~/paper_alleles/analysis/plots/figS1.pdf",figS1,units = "in",width = 8,height = 8,bg="white")
ggsave("~/paper_alleles/analysis/plots/figS1.png",figS1,units = "in",width = 8,height = 8,bg="white")
```

## FigS4c/d. SHM

```{r}
igenotyper_muts<-read_tsv(dir("~/paper_alleles/benchmarking/igenotyper/mixcr/mutations",
                              pattern=".tsv",
                              full.names = T),id = "filename")

igenotyper_muts <- igenotyper_muts |> 
  rename(shm=`nMutationsRate{FR1Begin(6):FR3End}`) |> 
  mutate(donor=str_remove_all(filename,".*\\/|_IGH.tsv|SRR.*?_"))

shmStats <- igenotyper_muts |> 
  filter(shm!="region_not_covered") |> 
  mutate(shm=as.numeric(shm)) |> 
  group_by(donor) |> 
  summarise(meanSHM=mean(shm))

g1<-comparisonStats %>%
    filter(
        totalClonotypes > 3000,
        matchStatus %in% c("falseNegative"),
        tool=="mixcr"
    ) %>%
    left_join(shmStats, by = "donor") %>%
    ggplot(aes(y = fraction, x = meanSHM)) +
    geom_point() +
    stat_cor(aes(label = paste(..r.label.., ..p.label.., sep = "~`,`~")), method = "pearson",
             label.x.npc = "middle", label.y.npc = "top") +
    theme_bw()+
   labs(y="Fraction of\nmissed alleles",x = "Mean SHM level in a sample")

g2<-comparisonStats %>%
    filter(
        totalClonotypes > 3000,
        matchStatus %in% c("falsePositive"),
        tool=="mixcr"
    ) %>%
    mutate(matchStatus = factor(matchStatus, levels = c("falseNegative", "match", "falsePositive"))) %>%
    left_join(shmStats, by = "donor") %>%
    ggplot(aes(y = N, x = meanSHM)) +
    geom_point() +
    stat_cor(aes(label = paste(..r.label.., ..p.label.., sep = "~`,`~")), method = "pearson",
             label.x.npc = "middle", label.y.npc = "top") +
    theme_bw()+
   labs(y="# potential false\npositive calls",x = "Mean SHM level in a sample")

```

## FigS4 A/B. UMI vs no-UMI

```{r}

facet_labels <- c(mixcr = "MiXCR", tigger = "TIgGER")

g3<-comparisonStats  %>% 
    left_join(protocol_donor,by="donor") |> 
    filter(
        totalClonotypes>3000,
        matchStatus %in% c("match"),
        tool=="mixcr") %>%   
    mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) |> 
    ggplot(aes(y=fraction,x=protocol))+
    geom_boxplot(fill="#198020")+
  scale_x_discrete(labels = c("5race" = "With UMIs", "noumi" = "Without UMIs")) +
    stat_compare_means(method = "kruskal", label.x = 1.35, label.y = 0.0)+
  theme_bw()+
  labs(y="Fraction of alleles\nsupported by long reads")

g4<-comparisonStats  %>% 
  left_join(protocol_donor,by="donor") |> 
  filter(
    totalClonotypes>3000,
    matchStatus %in% c("falsePositive"),
        tool=="mixcr") %>%   
  mutate(matchStatus=factor(matchStatus,levels=c("falseNegative","match","falsePositive"))) %>% 
  ggplot(aes(y=N,x=protocol))+
    geom_boxplot(fill="#198020")+
  scale_x_discrete(labels = c("5race" = "With UMIs", "noumi" = "Without UMIs")) +
    stat_compare_means(method = "kruskal", label.x = 1.5, label.y = 7.5)+
  theme_bw()+
 
  labs(y="# potential false\npositive calls")

```

## FigS4

```{r}


figS4<-plot_grid(
  plot_grid(g3,g4,ncol = 2,labels=c("A","B")),
  plot_grid(g1,g2,nrow = 2,labels=c("C","D")),
  rel_heights = c(0.7,1.2),nrow=2
  )
                  

ggsave(filename = "plots/figS4.pdf", plot= figS4,device = "pdf",
       width = 8,height =6,units = "in",bg="white")
ggsave(filename = "plots/figS4.png", plot= figS4,device = "png",
       width = 8,height = 8,units = "in",bg="white")
```

## Citations

```{r}
grateful::cite_packages(out.dir = ".",citation.style = "apa-6th-edition")
```

## FigS3. Time benchmarking plot

```{r}
donorsForBench<-comparisonStats |> filter(tool=='mixcr',depth=='full',totalClonotypes>3000) |> pull(donor) |> unique() 

names(bench_full)<-dir("~/paper_alleles/benchmarking/immcantation_fasta",
                  pattern="_db-pass.tsv",
                  full.names = T) |> 
  str_extract("S?[WC]-[0-9]{1,3}")

figS3<-bench_full[names(bench_full) %in% donorsForBench] |> 
    keep( ~ is_tibble(.x)) |> bind_rows(.id="donor") |> 
    mutate(total_time.tigger = as.numeric(as_bench_time(total_time))) |> 
  left_join(
    c(
      dir("~/paper_alleles/benchmarking/mixcr/5race/findAlleles", pattern = "time.log",full.names = T),
      dir("~/paper_alleles/benchmarking/mixcr/noumi/findAlleles", pattern = "time.log",full.names = T)
    ) |> 
      read_tsv(id='filename',col_names = c('empty','real','user',"sys","avmem","maxmem")) |> 
      mutate(donor=filename |>  str_extract("S?[WC]-[0-9]{1,3}"),
             total_time.mixcr=str_remove(real," real,") |> as.numeric()) |> 
      select(donor,total_time.mixcr),
    by="donor"
  ) |> 
  pivot_longer(cols=starts_with("total_time."),names_to = "tool") |> 
  mutate(tool=str_remove(tool,"total_time.")) |> 
  ggplot(aes(x=tool,y=value,fill=tool))+
  geom_boxplot()+
  geom_point()+
  theme_bw()+
  scale_fill_manual(values = mi_dark,labels=c("MiXCR","TIgGER"))+
  stat_compare_means(method="wilcox")+
  scale_x_discrete(labels=c(mixcr="MiXCR",  tigger="TIgGER"))+
  labs(y="Time (s)", x= "Tool")

ggsave(filename = "plots/figS3.pdf", plot= figS3,device = "pdf",
       width = 6,height =4,units = "in",bg="white")
ggsave(filename = "plots/figS3.png", plot= figS3,device = "png",
       width = 6,height = 4,units = "in",bg="white")
```
