B:BD[
6.8430] → [
6.8430:8750]
B:BD[
6.8750] → [
8.49503:62716]
checkIfExists: true)
bamIndex = bam.map { it -> it + ".bai" }
downloadGenome | indexGenome
indexGenome.out.index | view
addSNV(f, bam, bamIndex, downloadGenome.out, indexGenome.out.index, indexGenome.out.dict, indexGenome.out.fai)
}
#+end_src
****** DONE v1.3: Lancer le pipeline pour vérifie
r qu'on retrouve les variants
CLOSED: [2023-05-19 Fri 18:41] SCHEDULED: <2023-05-18 Thu>
****** TODO Corrigier position pour avoir une bonne VAF
SCHEDULED: <2023-05-19 Fri>
POS POS+1 VAF ALT
Attention, la base corrigée est à POS+1...
****** DONE Comparaison manuelle avec julia (VAF = 1...)
CLOSED: [2023-05-19 Fri 21:58]
552 found over 585
***** TODO del
***** TODO ins
*** TODO Julia : Bamscissors :bamscissors:
Modification de la séquence dans BAM.
*Pas de mise à jour de CIGAR*
On convertit en fastq et on lance le pipeline pour "corriger"
#+begin_src sh
cd /home/alex/code/bisonex/out/63003856/preprocessing/mapped
samtools view 63003856_S135.bam NC_000022.11 -o 63003856_S135_chr22.bam
cd /home/alex/recherche/bisonex/code/BamScissors.jl
cp ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135_chr22.bam .
samtools index 63003856_chr22.bam
#+end_src
Le script va modifier le bam, le trier et générer le fastq. !!!
Attention: ne pas oublier l'option -n !!!
#+begin_src sh
time julia --project=.. insertVariant.jl
scp 63003856_S135_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
**** TODO Phase 1 : chr22, VAF=1
***** DONE 1 SNV : ok !
CLOSED: [2023-05-20 Sat 19:35] SCHEDULED: <2023-05-20 Sat>
#+begin_src
make run READS="tests/bamscissors/corrected_{1,2}.fq.gz"
Puis on lance le pipeline sur correct_1
- [X] Variant visible dans IGV
- [X] Variant visible après alignement
- [X] Variant visible après appel de variant
***** TODO Tous les SNV du chromosome
SCHEDULED: <2023-05-20 Sat>
**** TODO PHase 2 : chr22, VAF variable: de nombreux reads sont perdus
Problème dansr le BAM car sans les insertion
***** DONE Filtrer les reads sans pair : idem
CLOSED: [2023-05-23 Tue 00:01]
FOUND:Found 16662 unpaired mates
at htsjdk.samtools.SAMUtils.processValidationError(SAMUtils.java:470)
at picard.sam.SamToFastq.doWork(SamToFastq.java:224)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
ex: chr22:42213078
On essaie
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
Ne rale pas...
SCHEDULED: <2023-05-22 Mon>
***** TODO Problème de la conversion BAM -> fastq ?
****** DONE Test bamtofastq sur le BAM originel: idem
CLOSED: [2023-05-23 Tue 01:16]
Dans ~/recherche/code/bisonex/Bamscissors.jl
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
On trie bien par nom de read
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
Conversion en fastq. Attention à bien prendre le *fichier trié* !
samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
On envoie sur le mésocentre
scp 63003856_chr22_init{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors
On démarre un job avec 24 coeurs pour aller vite
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_init1.fq.gz 63003856_chr22_init2.fq.gz | samtools sort -@24 -o output.bam -
Dans /home/alex/recherche/bisonex/code/BamScissors.jl/aligned
❯ scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam*
****** DONE Avec samtools fastq
CLOSED: [2023-05-23 Tue 01:16]
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -o 63003856_chr22_sorted.bam
scp 63003856_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
mesocentre
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam\* aligned/
****** DONE En rajoutant .fna dans le doserrie (+samtools fastq): idem
CLOSED: [2023-05-23 Tue 01:16]
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef* .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef* .
bwa mem -t 24 genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
****** DONE Test picard : idem
CLOSED: [2023-05-23 Tue 01:16]
Dans bisonex/code/BamScissors.jl
❯ picard SamToFastq -I 63003856_chr22_sorted.bam -F testpicard1.fastq -F2 testpicard2.fastq
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef testpicard1.fastq.gz testpicard2.fastq.gz | samtools sort -@24 -o testpicard.bam -
****** DONE Il ne manque pas des reads dans les fastq :
CLOSED: [2023-05-23 Tue 01:26]
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat aligned/output.bam
1306273 + 0 in total (QC-passed reads + QC-failed reads)
1296055 + 0 primary
0 + 0 secondary
10218 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
1304871 + 0 mapped (99.89% : N/A)
1294653 + 0 primary mapped (99.89% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
❯ echo $(zcat 63003856_chr22_init2.fq.gz | wc -l)/4 | bc
1296055
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0 took 2s
❯ echo $(zcat 63003856_chr22_init1.fq.gz | wc -l)/4 | bc
1296055
❯ samtools flagstat 63003856_chr22_sorted.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
****** DONE Avec sort -O BAM idem
CLOSED: [2023-05-23 Tue 01:21]
****** DONE Singleton ou non mapped ? nonVirtPosition
CLOSED: [2023-05-23 Tue 01:22]
❯ samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -0 both -s single.fq.gz -n 63003856_chr22_sorted.bam
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ wc -l 63003856_chr22_init* both single.fq.gz
403540 63003856_chr22_init1.fq.gz
404788 63003856_chr22_init2.fq.gz
0 both
0 single.fq.gz
808328 total
****** Problème d'aligner car les reads sont bien dans le .fastq ?
Ex:
zgrep "A00853:477:HMLWYDSX3:2:2624:2826:18630" 63003856_chr22_{1,2}.fq.gz
63003856_chr22_1.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
63003856_chr22_2.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
***** Refaire la manip sur bam chr22 non modifié + mail Alexis
[1]_samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
Les reads sont bien tous mappé
samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
[2] samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
flagstat idem
[3] samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
Nombre de reads ok (= la moitié) et les séquences ne sont pas identiques pour le premier read (= on n'a pas 2x la même chose)
echo $(zcat 63003856_chr22_1.fq.gz | wc -l)/4 | bc
1296055
echo $(zcat 63003856_chr22_2.fq.gz | wc -l)/4 | bc
1296055
[4]
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz -o wtf.bam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1752492 sequences (240000014 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (12, 702821, 0, 18)
[M::mem_pestat] analyzing insert size distribution for orientation FF...
[M::mem_pestat] (25, 50, 75) percentile: (54, 197, 268)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 696)
[M::mem_pestat] mean and std.dev: (163.92, 129.71)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 910)
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (129, 175, 233)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 441)
[M::mem_pestat] mean and std.dev: (185.56, 75.37)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 545)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (56, 192, 239)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 605)
[M::mem_pestat] mean and std.dev: (158.00, 110.30)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 788)
[M::mem_pestat] skip orientation FF
[M::mem_pestat] skip orientation RR
[M::process] read 839618 sequences (114952336 bp)...
[M::mem_process_seqs] Processed 1752492 reads in 375.714 CPU sec, 17.645 real sec
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 336379, 0, 5)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (128, 174, 232)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 440)
[M::mem_pestat] mean and std.dev: (184.73, 74.63)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 544)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 839618 reads in 183.039 CPU sec, 7.961 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 24 -o wtf.bam /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz
[main] Real time: 38.278 sec; CPU: 565.821 sec
Bon nombre de reads pourtant
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools sort -@24 -o wtf_sorted.bam wtf.sam
[bam_sort_core] merging from 0 files and 24 in-memory blocks...
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
Effectivement, ce n'est pas un problème d'IGV
❯ samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 19 TTtTTtTTTtTTtTtTttt kkk_kkkkFkFF_FkFkQk
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools mpileup aligned/wtf_sorted.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 5 TTtTT _FkFF
*** Divers
**** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
checkIfExists: true)
bamIndex = bam.map { it -> it + ".bai" }
downloadGenome | indexGenome
indexGenome.out.index | view
addSNV(f, bam, bamIndex, downloadGenome.out, indexGenome.out.index, indexGenome.out.dict, indexGenome.out.fai)
}
#+end_src
****** DONE v1.3: Lancer le pipeline pour vérifier qu'on retrouve les variants
CLOSED: [2023-05-19 Fri 18:41] SCHEDULED: <2023-05-18 Thu>
****** TODO Corrigier position pour avoir une bonne VAF
SCHEDULED: <2023-05-19 Fri>
POS POS+1 VAF ALT
Attention, la base corrigée est à POS+1...
****** DONE Comparaison manuelle avec julia (VAF = 1...)
CLOSED: [2023-05-19 Fri 21:58]
552 found over 585
***** TODO del
***** TODO ins
*** TODO Julia : Bamscissors :bamscissors:
Modification de la séquence dans BAM.
*Pas de mise à jour de CIGAR*
On convertit en fastq et on lance le pipeline pour "corriger"
#+begin_src sh
cd /home/alex/code/bisonex/out/63003856/preprocessing/mapped
samtools view 63003856_S135.bam NC_000022.11 -o 63003856_S135_chr22.bam
cd /home/alex/recherche/bisonex/code/BamScissors.jl
cp ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135_chr22.bam .
samtools index 63003856_chr22.bam
#+end_src
Le script va modifier le bam, le trier et générer le fastq. !!!
Attention: ne pas oublier l'option -n !!!
#+begin_src sh
time julia --project=.. insertVariant.jl
scp 63003856_S135_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
**** TODO Phase 1 : chr22, VAF=1
***** DONE 1 SNV : ok !
CLOSED: [2023-05-20 Sat 19:35] SCHEDULED: <2023-05-20 Sat>
#+begin_src
make run READS="tests/bamscissors/corrected_{1,2}.fq.gz"
Puis on lance le pipeline sur correct_1
- [X] Variant visible dans IGV
- [X] Variant visible après alignement
- [X] Variant visible après appel de variant
***** TODO Tous les SNV du chromosome
SCHEDULED: <2023-05-20 Sat>
**** DONE PHase 2 : chr22, VAF variable: de nombreux reads sont perdus -> ok sur un SNV en alignant sur chromosome 22
CLOSED: [2023-05-24 Wed 23:19]
Problème dansr le BAM car sans les insertion
***** DONE Filtrer les reads sans pair : idem
CLOSED: [2023-05-23 Tue 00:01]
FOUND:Found 16662 unpaired mates
at htsjdk.samtools.SAMUtils.processValidationError(SAMUtils.java:470)
at picard.sam.SamToFastq.doWork(SamToFastq.java:224)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
ex: chr22:42213078
On essaie
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
Ne rale pas...
SCHEDULED: <2023-05-22 Mon>
***** DONE Problème de la conversion BAM -> fastq ?
CLOSED: [2023-05-24 Wed 23:19]
****** DONE Test bamtofastq sur le BAM originel: idem
CLOSED: [2023-05-23 Tue 01:16]
Dans ~/recherche/code/bisonex/Bamscissors.jl
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
On trie bien par nom de read
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
Conversion en fastq. Attention à bien prendre le *fichier trié* !
samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
On envoie sur le mésocentre
scp 63003856_chr22_init{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors
On démarre un job avec 24 coeurs pour aller vite
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_init1.fq.gz 63003856_chr22_init2.fq.gz | samtools sort -@24 -o output.bam -
Dans /home/alex/recherche/bisonex/code/BamScissors.jl/aligned
❯ scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam*
****** DONE Avec samtools fastq
CLOSED: [2023-05-23 Tue 01:16]
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -o 63003856_chr22_sorted.bam
scp 63003856_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
mesocentre
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam\* aligned/
****** DONE En rajoutant .fna dans le doserrie (+samtools fastq): idem
CLOSED: [2023-05-23 Tue 01:16]
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef* .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef* .
bwa mem -t 24 genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
****** DONE Test picard : idem
CLOSED: [2023-05-23 Tue 01:16]
Dans bisonex/code/BamScissors.jl
❯ picard SamToFastq -I 63003856_chr22_sorted.bam -F testpicard1.fastq -F2 testpicard2.fastq
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef testpicard1.fastq.gz testpicard2.fastq.gz | samtools sort -@24 -o testpicard.bam -
****** DONE Il ne manque pas des reads dans les fastq :
CLOSED: [2023-05-23 Tue 01:26]
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat aligned/output.bam
1306273 + 0 in total (QC-passed reads + QC-failed reads)
1296055 + 0 primary
0 + 0 secondary
10218 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
1304871 + 0 mapped (99.89% : N/A)
1294653 + 0 primary mapped (99.89% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
❯ echo $(zcat 63003856_chr22_init2.fq.gz | wc -l)/4 | bc
1296055
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0 took 2s
❯ echo $(zcat 63003856_chr22_init1.fq.gz | wc -l)/4 | bc
1296055
❯ samtools flagstat 63003856_chr22_sorted.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
****** DONE Avec sort -O BAM idem
CLOSED: [2023-05-23 Tue 01:21]
****** DONE Singleton ou non mapped ? nonVirtPosition
CLOSED: [2023-05-23 Tue 01:22]
❯ samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -0 both -s single.fq.gz -n 63003856_chr22_sorted.bam
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ wc -l 63003856_chr22_init* both single.fq.gz
403540 63003856_chr22_init1.fq.gz
404788 63003856_chr22_init2.fq.gz
0 both
0 single.fq.gz
808328 total
****** Problème d'aligner car les reads sont bien dans le .fastq ?
Ex:
zgrep "A00853:477:HMLWYDSX3:2:2624:2826:18630" 63003856_chr22_{1,2}.fq.gz
63003856_chr22_1.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
63003856_chr22_2.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
***** DONE Refaire la manip sur bam chr22 non modifié + mail Alexis
CLOSED: [2023-05-23 Tue 23:27]
[1]_samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
Les reads sont bien tous mappé
samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
[2] samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
flagstat idem
[3] samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
Nombre de reads ok (= la moitié) et les séquences ne sont pas identiques pour le premier read (= on n'a pas 2x la même chose)
echo $(zcat 63003856_chr22_1.fq.gz | wc -l)/4 | bc
1296055
echo $(zcat 63003856_chr22_2.fq.gz | wc -l)/4 | bc
1296055
[4]
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz -o wtf.bam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1752492 sequences (240000014 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (12, 702821, 0, 18)
[M::mem_pestat] analyzing insert size distribution for orientation FF...
[M::mem_pestat] (25, 50, 75) percentile: (54, 197, 268)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 696)
[M::mem_pestat] mean and std.dev: (163.92, 129.71)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 910)
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (129, 175, 233)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 441)
[M::mem_pestat] mean and std.dev: (185.56, 75.37)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 545)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (56, 192, 239)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 605)
[M::mem_pestat] mean and std.dev: (158.00, 110.30)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 788)
[M::mem_pestat] skip orientation FF
[M::mem_pestat] skip orientation RR
[M::process] read 839618 sequences (114952336 bp)...
[M::mem_process_seqs] Processed 1752492 reads in 375.714 CPU sec, 17.645 real sec
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 336379, 0, 5)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (128, 174, 232)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 440)
[M::mem_pestat] mean and std.dev: (184.73, 74.63)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 544)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 839618 reads in 183.039 CPU sec, 7.961 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 24 -o wtf.bam /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz
[main] Real time: 38.278 sec; CPU: 565.821 sec
Bon nombre de reads pourtant
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools sort -@24 -o wtf_sorted.bam wtf.sam
[bam_sort_core] merging from 0 files and 24 in-memory blocks...
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
Effectivement, ce n'est pas un problème d'IGV
❯ samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 19 TTtTTtTTTtTTtTtTttt kkk_kkkkFkFF_FkFkQk
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools mpileup aligned/wtf_sorted.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 5 TTtTT _FkFF
***** DONE Regarder où ont été aligné les reads (nouveau run)
CLOSED: [2023-05-24 Wed 23:18]
****** DONE Préparation
CLOSED: [2023-05-24 Wed 21:59]
On relance le pipeline pour avoir un BAM propre
On garde les reads non mappé à partsir de la sortie d'applybqsr
#+begin_src sh
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios --input="/Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R{1,2}_001.fastq.gz" --outdir=out -bg
cd out/63003856_S135_R/preprocessing/applybqsr/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
make run BG= READS="out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz"
cd out/63003856_chr22/preprocessing/mapped/
samtools index 63003856_chr22.bam
samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
#+end_src
On récupère les 2 bam dans
#+begin_src
cd /home/alex/recherche/bisonex/code/BamScissors.jl/
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_chr22/preprocessing/mapped data
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing/applybqsr/ data/init/
#+end_src
****** Vérification que le reads est ailleurs
On cherche un read manquant dans le second alignement
#+begin_src sh
samtools view data/init/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
#+begin_src sh
samtools view data/mapped/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
Effectivement, on aligne sur une zonne supprimée !
***** DONE Corriger la qualité: non
CLOSED: [2023-05-24 Wed 22:19]
****** DONE Comparaison avec le fastq de référénce : qualité !!
CLOSED: [2023-05-24 Wed 22:17]
#+begin_src sh
cd /Work/Users/apraga/bisonex/work/6e/8548fc90263830bf677f36585f11dc
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
@A00853:477:HMLWYDSX3:1:1413:4390:28573
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
+
ADEEB@?CBBCCBDCBDCCCFBD;EEBEBBABCEC:EEBEAADCFCDFBFECCFCFFFFCCGFEDGFBBEEDCAECCEEBCECBDBCEEDCCBCBFAECCFEACAEAEBCCDCBCBFB:;CAEDCAEDBEEEEDC?<ECFBCDBCCCEDAA
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" /Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R1_001.fastq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573 1:N:0:ATTCCACACA+TAGGCGATTG
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
****** DONE Regarder la qualité après bwa mem vs applybqsr: différente
CLOSED: [2023-05-24 Wed 22:18]
Sur le mésocentre, dans /Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing
$ samtools view mapped/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
samtools view applybqsr/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
****** DONE Réaligner à partir de la sortie de bwa mem
CLOSED: [2023-05-24 Wed 22:32]
#+begin_src sh
cd out/63003856_S135_R/preprocessing/mapped/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
#+end_src
ON vérifie la qualité
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573
: AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
#+begin_src
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios -resume --input="out/63003856_S135_R/preprocessing/mapped/63003856_chr22_{1,2}.fq.gz" --outdir=out/63003856_chr22-from-mapped
#+end_src
Puis ::
#+begin_src
cd /Work/Users/apraga/bisonex/out/63003856_chr22-from-mapped/63003856_chr22/preprocessing/mapped
samtools view 63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
***** DONE Aligner sur génome de référence limité au chromosome 22
CLOSED: [2023-05-24 Wed 23:18]
****** KILL Test données non modifiées
CLOSED: [2023-05-24 Wed 23:18]
/Work/Users/apraga/bisonex/tests/bamscissors
#+begin_src
cd /Work/Groups/bisonex/data/genome/GRCh38.p13/
mkdir chr22/
samtools faidx genomeRef.fna NC_000022.11 > chr22/chr22.fna
cd chr22
samtools faidx chr22.fna
bwa index chr22.fna
#+end_src
#+begin_src
cd /Work/Users/apraga/bisonex/tests/bamscissors
ln -s ../../out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz .
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz -o smallref.sam
#+end_src
****** DONE Test données modifiées: ok
CLOSED: [2023-05-24 Wed 23:18]
Données dans data/init
#+begin_src sh
time julia insertVariant.jl
rsync -avz data/init/*.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
#+begin_src
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz | samtools sort -@24 - -o smallref.bam
#+end_src
#+begin_src
rsync -avz meso:/Work/Users/apraga/bisonex/tests/bamscissors/smallref.bam mapped/
#+end_src
**** TODO Phase 3 : insertion SNV chromosome 22 sur BAM complet + test haplotypecaller
*** Divers
**** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]