Embedded Files
Environment setting @R:
Environment setting @R:
library(dada2);library(Biostrings); library(ShortRead); # From Bioconductor
library(ggplot2); library(reshape2); library(gridExtra) # installed from CRAN
DADA2, Biostrings, ShortRead
DADA2, Biostrings, ShortRead
Install
Install
ref. DADA2 & Bioconductor
Install by Binary from Bioconductor
@R$ if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")@R$ BiocManager::install("dada2", version = "3.21") # for R ver. 4.5#@R$ BiocManager::install("dada2", version = "3.19") # for R ver. 4.4@R$ BiocManager::install("Biostrings")@R$ BiocManager::install("ShortRead")DATA：
DATA：
PacBio-16S
PacBio-16S
Technology = PacBio SeqWell2
Amplicon = 16S rRNA gene long read with Primers（27F/1492R）
- F27:5'-AGRGTTYGATYATGGCTCAG-3'
- R1492:5'-RGTTACCTTGTTACGACTT-3'(R1492-Revcom: AAGTCGTAACAAGGTAACY)
生物技研使用プライマー：
- 16S_Forward：/5Phos/GCATC-[PacBio barcode sequences]-AGRGTTYGATYMTGGCTCAG (RC-For: CTGAGCCAKMATCMAACYCT)
- 16S_Reverse：/5Phos/GCATC-[PacBio barcode sequences]-RGYTACCTTGTTACGACTT (RC-Rev: AAGTCGTAACAAGGTARCY)
input file = [Unziped.fastq]
PacBio-18S
PacBio-18S
ref. benjjneb.github
まず、.fastqのままプライマーを除く必要がある！
まず、.fastqのままプライマーを除く必要がある！
いずれの有力なtoolもアダプター削除用なので縮重塩基には弱い。
Primer部分を認識しているようだが、プライマー部分をカットしない
残ったプライマー領域を対象に何通りも繰り返す必要がある。
▶︎ cutadapt ... 縮重塩基には非対応
NNN+Primerも改善せず。。。
PrimerのReverseComplementを考慮して与える必要あり
% cutadapt -g ^AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -o mao7_trm_max16_cutp1.fastq.gz mao7_trmd_max1600.fastq.gz
▶︎ seqkit amplicon .... 縮重塩基には非対応だが「ミスマッチ許容無し」設定でプライマー内ミスマッチを除外するためプライマーの長さが統一される。そのため、fastpにより塩基数を指定して前後の配列を除去することでプライマー領域を取り除くことができる。
--max-mismatchと縮重プライマーを同時に使用すると-mオプションは無視されるらしい (Ushio's)
←そうでもなさそう。-mの数字により結果が異なる
NNN+Primerも改善せず。。。
PrimerのReverseComplementを考慮する必要なし（F, Rともに5'→3'）
-m 4で良さそう
プライマーミスマッチを除去 & 断片の方向を揃える
% seqkit amplicon
-F AGRGTTYGATYMTGGCTCAG
-R RGYTACCTTGTTACGACTT # -m 4 ⬅︎ max-mismatchオプション無し=ミスマッチ許容しない
mao7_trmd_max1600.fastq.gz
> mao7_trm_m16_rmp_m4v10.fastq
fastpによりプライマーの塩基数を両端から削除する
% fastp -i [FILE: mao7.fastq] -o [FILE: mao7_rmfr.fastq] -f 20 -t 19
▶︎ clsplitseq (in claident) ... 縮重塩基には非対応
$ clsplitseq --runname=MAO7 --indexname=M --primerfile=/home/tak-gk/local/share/primers/Fprimer_BACFULL.fas --tagname=M --reverseprimerfile=/home/tak-gk/local/share/primers/Rprimer_BACFULL2.fas --truncateN=ENABLE --append --numthreads=4 ~/wd/WholeGenome/MAO7/Pac16S/mao7_trmd_max1600.fastq.gz ./1_rmprimers_Rev2Primer-trimming Log
Primer-trimming Log
▶︎ cutadapt
▶︎ cutadapt
... 縮重塩基には非対応
NNN+Primerも改善せず。。。
1st
cutadapt
~/wd/WholeGenome/MAO7/Pac16S$
% cutadapt
-g ^AGRGTTYGATYMTGGCTCAG
-a AAGTCGTAACAAGGTARCY
-o mao7_trm_max16_cutp1.fastq.gz
mao7_trmd_max1600.fastq.gz
2nd_Reverse case
~/wd/WholeGenome/MAO7/Pac16S$
% cutadapt
-g ^RGYTACCTTGTTACGACTT
-a CTGAGCCAKRATCRAACYCT
-o mao7_trm_max16_cutp_2nd.fastq.gz
mao7_trm_max16_cutp.fastq.gz
$ cutadapt -g AGRGTTYGATYMTGGCTCAG [F] -a AAGTCGTAACAAGGTARCY [Rrc]-e 0.1 -o mao7_trmd_max16_cutp1.fastq ➡︎ ../mao7_trmd_max1600.fastq.gz ⬅︎ 98,308 sequences (1,479–1,591 bp)
$ cutadapt -g CTGAGCCAKMATCMAACYCT [Frc] -a RGYTACCTTGTTACGACTT [R] -e 0.1 -o mao7_trmd_max16_cutp2.fastq ./mao7_trmd_max16_cutp1.fastq
$ cutadapt -g CTGAGCCAKMATCMAACYCT [Frc] -a RGYTACCTTGTTACGACTT [R] -e 0.1 -m 100 -o mao7_trmd_max16_cutp2.fastq ./mao7_trmd_max16_cutp1.fastq
$ cutadapt -g AGRGTTYGATYMTGGCTCAG [F]-e 0.1 -o mao7_trmd_max16_cutp1b.fastq ../mao7_trmd_max1600.fastq.gz
$ cutadapt -g AGRGTTYGATYMTGGCTCAG [F] -a AAGTCGTAACAAGGT [Rrc] -e 0.1 -o mao7_trmd_max16_cutp1b.fastq ../mao7_trmd_max1600.fastq.gz
$ cutadapt -g AGRGTTYGATYMTGGCTCAG [F] -a AAGTCGTAACAAGGT [Rrc]-e 0.2 -o mao7_trmd_max16_cutp1b.fastq ../mao7_trmd_max1600.fastq.gz
$ cutadapt -g ^.*AGRGTTYGATYMTGGCTCAG [F] -a AAGTCGTAACAAGGT.*$ [Rrc] -e 0.2 -o mao7_trmd_max16_cutp1c.fastq ../mao7_trmd_max1600.fastq.gz
~/wd/WholeGenome/MAO7/Pac16S/wd3/1_rmprimers_Rev2$ cutadapt -g RGYTACCTTGTTACGACTT [R] -a CTGAGCCAKMATCMAACYCT [Frc] -e 0.2 -m 100 -o mao7_trmd_max16_cutp2b.fastq ./mao7_trmd_max16_cutp1b.fastq
▶︎ seqkit amplicon
▶︎ seqkit amplicon
... --max-mismatchと縮重プライマーを同時に使用すると-mオプションは無視されるらしい (Ushio's)
$ seqkit amplicon
-F AGRGTTYGATYMTGGCTCAG
-R RGYTACCTTGTTACGACTT
-m 4
mao7_trmd_max1600.fastq.gz
> mao7_trm_m16_rmp_m4v10.fastq
▶︎ clsplitseq
▶︎ clsplitseq
$ clsplitseq --runname=MAO7 --indexname=M --primerfile=/home/tak-gk/local/share/primers/Fprimer_BACFULL.fas --tagname=M --reverseprimerfile=/home/tak-gk/local/share/primers/Rprimer_BACFULL2.fas --truncateN=ENABLE --append --numthreads=4 ~/wd/WholeGenome/MAO7/Pac16S/mao7_trmd_max1600.fastq.gz ./1_rmprimers_Rev2
R環境設定
R@ library(dada2); library(Biostrings); library(ShortRead);
library(ggplot2); library(reshape2); library(gridExtra); library(phyloseq) # installed from CRAN
Dereplicate:
R@ drp2 <- derepFastq('../working/', verbose=T)
- Dereplicating sequence entries in Fastq file: ../working//MAO6_LR16S_keep75.fastq
- Encountered 36055 unique sequences from 61498 total sequences read.
Learning errors
R@ err2 <- learnErrors(drp2, errorEstimationFunction=PacBioErrfun, BAND_SIZE=32, multithread=T)
- 93761581 total bases in 61498 reads from 1 samples will be used for learning the error rates.
R@ saveRDS(err2, "./MAO6_err2.rds")
R@ plotErrors(err2)
Denoising
R@ dd2 <- dada(drp2, err=err2, BAND_SIZE=32, multithread=T)
- Sample 1 - 61498 reads in 36055 unique sequences.
R@ saveRDS(dd2, "MAO6_dd2.rds")
Wrightout sequence table
R@ st2 <- makeSequenceTable(dd2); dim(st2)
-----
Fasta形式へ変換
% cat dada_keep75v2.txt | sed "s/^/>ASV_/" | perl -pe 's/\s/\n/g' > MAO6kp75_dada2.fas
アラインメント（F/Rの方向を一方へ統一: --adjustdirection）
% mafft --thread 4 --adjustdirection MAO6kp75_dada2.fas > MAO6kp75_dada2_ali.fas
========== LOGs ==========
240831 Jn10c07-18S-Pacbio
> drp2 <- derepFastq("./Jn10c07L18S_keep90_min900_max2700.fq.gz", verbose=T)Dereplicating sequence entries in Fastq file: ./Jn10c07L18S_keep90_min900_max2700.fq.gzEncountered 43951 unique sequences from 91626 total sequences read.
> err2 <- learnErrors(drp2, errorEstimationFunction = PacBioErrfun, BAND_SIZE=32,multithread=T)239042957 total bases in 91626 reads from 1 samples will be used for learning the error rates.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.
> dd2 <- dada(drp2, err=err2,BAND_SIZE=32, multithread=T)Sample 1 - 91626 reads in 43951 unique sequences.> saveRDS(dd2,"./res/Jn10c07SSU_dd2.rds")
Sample 1 - 91626 reads in 43951 unique sequences.
PacBio-18S
PacBio-18S
Jn10c07 & Jn10c08
Jn10c07 & Jn10c08
QC
QC
fastp: 統合QCツール（2 colors対応）GitHub, Kazmax's web-fastp, 本家manual
QCレポート：
% fastp -i [FILE1] -3 -o [FILE1_res]
-h [report .html] -j [report.json] -q 15 -n 10 -w 4
-q=平均Q-val-n=１readあたりのNの数。設定値を超えたreadとpair-readを除去-t=read1のtailから除去する配列数[default=0]-T=指定無しでread2に対して-tに準ずる-f=read1のfrontから除去する配列数[default=0]-F=指定無しでread2に対して-tに準ずる-l=最低リード長。設定以下のreadを除去-w=thread数Seqkit stats
$ seqkit stats Jn10c07L18S.fq.gz
file format type num_seqs sum_len min_len avg_len max_len
Jn10c07L18S.fq.gz FASTQ DNA 103,496 265,603,042 51 2,566.3 5,272
Filtlong
% filtlong --keep_percent 90 --min_length 2000 hiyc8ont_raw.fq.gz | pigz > h8trmd.fq.gz
鎖長分布
$ perl ~/local/bin/fastq_avelength_morethan1k.perl Jn10c07L18S.fastq$ perl ~/local/bin/fastq_avelength_morethan1k.perl Jn10c08L18S.fastq
MAO7-16S-PacBio
MAO7-16S-PacBio
cutadapt_log
cutadapt_log
#==============================================250616 Again : because Previous try could not remove all primers1st round : Forward & RC-Reverse
$ cutadapt -g AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -e 0.1 -o mao7_trmd_max16_cutp1.fastq ../mao7_trmd_max1600.fastq.gz
=== Summary ===
Total reads processed: 98,308Reads with adapters: 62,684 (63.8%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 149,849,626 bpTotal written (filtered): 148,595,746 bp (99.2%)
=== Adapter 1 ===
Sequence: AGRGTTYGATYMTGGCTCAG; Type: regular 5'; Length: 20; Trimmed: 54799 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1; 20 bp: 2
Overview of removed sequenceslength count expect max.err error counts3 1 1536.1 0 118 1 0.0 1 119 335 0.0 1 23 31220 47899 0.0 2 47800 90 921 6167 0.0 2 12 615522 386 0.0 2 5 43 33823 1 0.0 2 0 125 1 0.0 2 148 2 0.0 2 251 1 0.0 2 154 1 0.0 2 155 1 0.0 2 156 1 0.0 2 160 1 0.0 2 161 1 0.0 2 1
=== Adapter 2 ===
Sequence: AAGTCGTAACAAGGTARCY; Type: regular 3'; Length: 19; Trimmed: 7885 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 0.0% C: 0.0% G: 100.0% T: 0.0% none/other: 0.0%WARNING: The adapter is preceded by 'G' extremely often. The provided adapter sequence could be incomplete at its 5' end. Ignore this warning when trimming primers.
Overview of removed sequenceslength count expect max.err error counts18 49 0.0 1 42 719 6630 0.0 1 6623 720 1169 0.0 1 979 19021 35 0.0 1 1 3422 1 0.0 1 145 1 0.0 1 1
# VVVVV !! Same logs as the first analysis (250603) VVVVVV
2nd round : Change primer order: RC-Forwerd & Reverse$ cutadapt -g CTGAGCCAKMATCMAACYCT -a RGYTACCTTGTTACGACTT -e 0.1 -m 100 -o mao7_trmd_max16_cutp2.fastq ./mao7_trmd_max16_cutp1.fastq
=== Summary ===
Total reads processed: 98,308Reads with adapters: 47,341 (48.2%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 148,595,746 bpTotal written (filtered): 94,246,003 bp (63.4%)
=== Adapter 1 ===
Sequence: CTGAGCCAKMATCMAACYCT; Type: regular 5'; Length: 20; Trimmed: 23573 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1; 20 bp: 2
Overview of removed sequenceslength count expect max.err error counts1500 2 0.0 2 21501 6 0.0 2 5 11502 4 0.0 2 41503 7 0.0 2 4 31504 15 0.0 2 9 5 11505 15 0.0 2 13 21506 30 0.0 2 24 61507 37 0.0 2 31 5 11508 59 0.0 2 47 121509 96 0.0 2 74 21 11510 140 0.0 2 111 28 11511 220 0.0 2 154 63 31512 429 0.0 2 344 81 41513 1464 0.0 2 1325 134 51514 1785 0.0 2 1499 275 111515 1100 0.0 2 733 343 241516 768 0.0 2 469 249 501517 557 0.0 2 333 178 461518 437 0.0 2 268 136 331519 373 0.0 2 228 121 241520 333 0.0 2 199 107 271521 361 0.0 2 240 106 151522 349 0.0 2 234 104 111523 380 0.0 2 264 107 91524 427 0.0 2 296 115 161525 417 0.0 2 287 110 201526 443 0.0 2 308 125 101527 528 0.0 2 375 141 121528 1006 0.0 2 799 186 211529 3233 0.0 2 2956 254 231530 3018 0.0 2 2298 677 431531 2072 0.0 2 1417 569 861532 1470 0.0 2 945 457 681533 994 0.0 2 580 349 651534 514 0.0 2 230 228 561535 260 0.0 2 115 96 491536 104 0.0 2 45 42 171537 57 0.0 2 18 27 121538 22 0.0 2 5 10 71539 14 0.0 2 8 4 21541 3 0.0 2 0 2 11542 5 0.0 2 2 2 11543 2 0.0 2 21544 1 0.0 2 0 11547 1 0.0 2 11548 1 0.0 2 11549 1 0.0 2 0 11552 2 0.0 2 1 11554 2 0.0 2 1 0 11557 1 0.0 2 0 11558 2 0.0 2 21559 1 0.0 2 0 0 11560 2 0.0 2 0 1 11568 1 0.0 2 0 11574 1 0.0 2 0 0 11580 1 0.0 2 0 1
=== Adapter 2 ===
Sequence: RGYTACCTTGTTACGACTT; Type: regular 3'; Length: 19; Trimmed: 23768 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 16.8% C: 0.0% G: 0.3% T: 38.6% none/other: 44.3%
Overview of removed sequenceslength count expect max.err error counts3 587 1536.1 0 5874 11128 384.0 0 111288 1 1.5 0 11500 3 0.0 1 31501 5 0.0 1 51502 3 0.0 1 31503 2 0.0 1 21504 11 0.0 1 111505 15 0.0 1 151506 15 0.0 1 14 11507 21 0.0 1 211508 36 0.0 1 34 21509 57 0.0 1 51 61510 81 0.0 1 78 31511 148 0.0 1 143 51512 280 0.0 1 267 131513 950 0.0 1 936 141514 1003 0.0 1 987 161515 442 0.0 1 434 81516 334 0.0 1 328 61517 247 0.0 1 243 41518 185 0.0 1 181 41519 174 0.0 1 170 41520 152 0.0 1 146 61521 119 0.0 1 1191522 166 0.0 1 165 11523 179 0.0 1 174 51524 214 0.0 1 210 41525 198 0.0 1 194 41526 230 0.0 1 221 91527 280 0.0 1 274 61528 691 0.0 1 675 161529 2202 0.0 1 2182 201530 1592 0.0 1 1566 261531 832 0.0 1 813 191532 611 0.0 1 598 131533 307 0.0 1 299 81534 137 0.0 1 134 31535 81 0.0 1 79 21536 21 0.0 1 19 21537 12 0.0 1 11 11538 8 0.0 1 7 11539 2 0.0 1 21541 1 0.0 1 11547 1 0.0 1 11550 1 0.0 1 0 11558 1 0.0 1 11559 1 0.0 1 11590 1 0.0 1 1
#==============================================tak-gk:~/wd/WholeGenome/MAO7/Pac16S$ cutadapt -g ^AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -o mao7_trm_max16_cutp.fastq.gz mao7_trmd_max1600.fastq.gz This is cutadapt 4.9 with Python 3.10.14Command line parameters: -g ^AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -o mao7_trm_max16_cutp.fastq.gz mao7_trmd_max1600.fastq.gzProcessing single-end reads on 1 core ...Done 00:00:03 98,308 reads @ 34.7 µs/read; 1.73 M reads/minuteFinished in 3.416 s (34.752 µs/read; 1.73 M reads/minute).
=== Summary ===
Total reads processed: 98,308Reads with adapters: 62,683 (63.8%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 149,849,626 bpTotal written (filtered): 148,606,661 bp (99.2%)
=== Adapter 1 ===
Sequence: AGRGTTYGATYMTGGCTCAG; Type: anchored 5'; Length: 20; Trimmed: 49496 times
No. of allowed errors: 2
Overview of removed sequenceslength count expect max.err error counts19 373 0.0 1 0 37320 47876 0.0 2 47800 75 121 1246 0.0 2 0 124622 1 0.0 2 0 0 1
=== Adapter 2 ===
Sequence: AAGTCGTAACAAGGTARCY; Type: regular 3'; Length: 19; Trimmed: 13187 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 0.0% C: 0.0% G: 100.0% T: 0.0% none/other: 0.0%WARNING: The adapter is preceded by 'G' extremely often. The provided adapter sequence could be incomplete at its 5' end. Ignore this warning when trimming primers.
Overview of removed sequenceslength count expect max.err error counts17 1 0.0 1 0 118 335 0.0 1 307 2819 11012 0.0 1 10991 2120 1803 0.0 1 1650 15321 31 0.0 1 4 2722 1 0.0 1 123 1 0.0 1 140 1 0.0 1 145 1 0.0 1 154 1 0.0 1 1
WARNING: One or more of your adapter sequences may be incomplete. Please see the detailed output above.
$ cutadapt -g AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -e 0.1 -o mao7_trmd_max16_cutp1.fastq ../mao7_trmd_max1600.fastq.gz
- -a (Reverse primer must be RC seq)
=== Summary ===
Total reads processed: 98,308Reads with adapters: 62,684 (63.8%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 149,849,626 bpTotal written (filtered): 148,595,746 bp (99.2%)
=== Adapter 1 ===
Sequence: AGRGTTYGATYMTGGCTCAG; Type: regular 5'; Length: 20; Trimmed: 54799 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1; 20 bp: 2
Overview of removed sequenceslength count expect max.err error counts3 1 1536.1 0 118 1 0.0 1 119 335 0.0 1 23 31220 47899 0.0 2 47800 90 921 6167 0.0 2 12 615522 386 0.0 2 5 43 33823 1 0.0 2 0 125 1 0.0 2 148 2 0.0 2 251 1 0.0 2 154 1 0.0 2 155 1 0.0 2 156 1 0.0 2 160 1 0.0 2 161 1 0.0 2 1
=== Adapter 2 ===
Sequence: AAGTCGTAACAAGGTARCY; Type: regular 3'; Length: 19; Trimmed: 7885 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 0.0% C: 0.0% G: 100.0% T: 0.0% none/other: 0.0%WARNING: The adapter is preceded by 'G' extremely often. The provided adapter sequence could be incomplete at its 5' end. Ignore this warning when trimming primers.
Overview of removed sequenceslength count expect max.err error counts18 49 0.0 1 42 719 6630 0.0 1 6623 720 1169 0.0 1 979 19021 35 0.0 1 1 3422 1 0.0 1 145 1 0.0 1 1
# VVVVV !! Same logs as the first analysis (250603) VVVVVV
2nd round : Change primer order: RC-Forwerd & Reverse$ cutadapt -g CTGAGCCAKMATCMAACYCT -a RGYTACCTTGTTACGACTT -e 0.1 -m 100 -o mao7_trmd_max16_cutp2.fastq ./mao7_trmd_max16_cutp1.fastq
=== Summary ===
Total reads processed: 98,308Reads with adapters: 47,341 (48.2%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 148,595,746 bpTotal written (filtered): 94,246,003 bp (63.4%)
=== Adapter 1 ===
Sequence: CTGAGCCAKMATCMAACYCT; Type: regular 5'; Length: 20; Trimmed: 23573 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1; 20 bp: 2
Overview of removed sequenceslength count expect max.err error counts1500 2 0.0 2 21501 6 0.0 2 5 11502 4 0.0 2 41503 7 0.0 2 4 31504 15 0.0 2 9 5 11505 15 0.0 2 13 21506 30 0.0 2 24 61507 37 0.0 2 31 5 11508 59 0.0 2 47 121509 96 0.0 2 74 21 11510 140 0.0 2 111 28 11511 220 0.0 2 154 63 31512 429 0.0 2 344 81 41513 1464 0.0 2 1325 134 51514 1785 0.0 2 1499 275 111515 1100 0.0 2 733 343 241516 768 0.0 2 469 249 501517 557 0.0 2 333 178 461518 437 0.0 2 268 136 331519 373 0.0 2 228 121 241520 333 0.0 2 199 107 271521 361 0.0 2 240 106 151522 349 0.0 2 234 104 111523 380 0.0 2 264 107 91524 427 0.0 2 296 115 161525 417 0.0 2 287 110 201526 443 0.0 2 308 125 101527 528 0.0 2 375 141 121528 1006 0.0 2 799 186 211529 3233 0.0 2 2956 254 231530 3018 0.0 2 2298 677 431531 2072 0.0 2 1417 569 861532 1470 0.0 2 945 457 681533 994 0.0 2 580 349 651534 514 0.0 2 230 228 561535 260 0.0 2 115 96 491536 104 0.0 2 45 42 171537 57 0.0 2 18 27 121538 22 0.0 2 5 10 71539 14 0.0 2 8 4 21541 3 0.0 2 0 2 11542 5 0.0 2 2 2 11543 2 0.0 2 21544 1 0.0 2 0 11547 1 0.0 2 11548 1 0.0 2 11549 1 0.0 2 0 11552 2 0.0 2 1 11554 2 0.0 2 1 0 11557 1 0.0 2 0 11558 2 0.0 2 21559 1 0.0 2 0 0 11560 2 0.0 2 0 1 11568 1 0.0 2 0 11574 1 0.0 2 0 0 11580 1 0.0 2 0 1
=== Adapter 2 ===
Sequence: RGYTACCTTGTTACGACTT; Type: regular 3'; Length: 19; Trimmed: 23768 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 16.8% C: 0.0% G: 0.3% T: 38.6% none/other: 44.3%
Overview of removed sequenceslength count expect max.err error counts3 587 1536.1 0 5874 11128 384.0 0 111288 1 1.5 0 11500 3 0.0 1 31501 5 0.0 1 51502 3 0.0 1 31503 2 0.0 1 21504 11 0.0 1 111505 15 0.0 1 151506 15 0.0 1 14 11507 21 0.0 1 211508 36 0.0 1 34 21509 57 0.0 1 51 61510 81 0.0 1 78 31511 148 0.0 1 143 51512 280 0.0 1 267 131513 950 0.0 1 936 141514 1003 0.0 1 987 161515 442 0.0 1 434 81516 334 0.0 1 328 61517 247 0.0 1 243 41518 185 0.0 1 181 41519 174 0.0 1 170 41520 152 0.0 1 146 61521 119 0.0 1 1191522 166 0.0 1 165 11523 179 0.0 1 174 51524 214 0.0 1 210 41525 198 0.0 1 194 41526 230 0.0 1 221 91527 280 0.0 1 274 61528 691 0.0 1 675 161529 2202 0.0 1 2182 201530 1592 0.0 1 1566 261531 832 0.0 1 813 191532 611 0.0 1 598 131533 307 0.0 1 299 81534 137 0.0 1 134 31535 81 0.0 1 79 21536 21 0.0 1 19 21537 12 0.0 1 11 11538 8 0.0 1 7 11539 2 0.0 1 21541 1 0.0 1 11547 1 0.0 1 11550 1 0.0 1 0 11558 1 0.0 1 11559 1 0.0 1 11590 1 0.0 1 1
#==============================================tak-gk:~/wd/WholeGenome/MAO7/Pac16S$ cutadapt -g ^AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -o mao7_trm_max16_cutp.fastq.gz mao7_trmd_max1600.fastq.gz This is cutadapt 4.9 with Python 3.10.14Command line parameters: -g ^AGRGTTYGATYMTGGCTCAG -a AAGTCGTAACAAGGTARCY -o mao7_trm_max16_cutp.fastq.gz mao7_trmd_max1600.fastq.gzProcessing single-end reads on 1 core ...Done 00:00:03 98,308 reads @ 34.7 µs/read; 1.73 M reads/minuteFinished in 3.416 s (34.752 µs/read; 1.73 M reads/minute).
=== Summary ===
Total reads processed: 98,308Reads with adapters: 62,683 (63.8%)Reads written (passing filters): 98,308 (100.0%)
Total basepairs processed: 149,849,626 bpTotal written (filtered): 148,606,661 bp (99.2%)
=== Adapter 1 ===
Sequence: AGRGTTYGATYMTGGCTCAG; Type: anchored 5'; Length: 20; Trimmed: 49496 times
No. of allowed errors: 2
Overview of removed sequenceslength count expect max.err error counts19 373 0.0 1 0 37320 47876 0.0 2 47800 75 121 1246 0.0 2 0 124622 1 0.0 2 0 0 1
=== Adapter 2 ===
Sequence: AAGTCGTAACAAGGTARCY; Type: regular 3'; Length: 19; Trimmed: 13187 times
Minimum overlap: 3No. of allowed errors:1-9 bp: 0; 10-19 bp: 1
Bases preceding removed adapters: A: 0.0% C: 0.0% G: 100.0% T: 0.0% none/other: 0.0%WARNING: The adapter is preceded by 'G' extremely often. The provided adapter sequence could be incomplete at its 5' end. Ignore this warning when trimming primers.
Overview of removed sequenceslength count expect max.err error counts17 1 0.0 1 0 118 335 0.0 1 307 2819 11012 0.0 1 10991 2120 1803 0.0 1 1650 15321 31 0.0 1 4 2722 1 0.0 1 123 1 0.0 1 140 1 0.0 1 145 1 0.0 1 154 1 0.0 1 1
WARNING: One or more of your adapter sequences may be incomplete. Please see the detailed output above.
1st
2nd remove primer
3rd remove primer accurately
4th dada2 ver 1.36
2025/6/3
Torimming
~/wd/WholeGenome/MAO7/Pac16S$ filtlong --keep_percent 90 --min_length 1500 MAO7-16S.fastq.gz | pigz > mao7_trmd.fq.gz
~/wd/WholeGenome/MAO7/Pac16S$ filtlong --keep_percent 90 --min_length 1500 --max_length 1600 MAO7-16S.fastq.gz | pigz > mao7_trmd_max1600.fq.gz
@R% drep2 <- derepFastq('./mao7_trmd_max1600.fastq.gz', verbose=T)Dereplicating sequence entries in Fastq file: ./mao7_trmd_max1600.fastq.gzEncountered 87128 unique sequences from 98308 total sequences read.@R% err2 <- learnErrors(drep2, errorEstimationFunction=PacBioErrfun, BAND_SIZE=32, multithread=T)149849626 total bases in 98308 reads from 1 samples will be used for learning the error rates.@R% saveRDS(err2,"./mao7_err2.rds")@R% plotErrors(err2)@R% dd2 <- dada(drep2,err=err2, BAND_SIZE=32,multithread=T)@R% saveRDS(dd2,"MAO7_dd2.rds")R@ st2 <- makeSequenceTable(dd2); dim(st2)
#-- rds to .fasta@R> seqtab <- readRDS("MAO7_dd2.rds")> asv_seqs <- colnames(seqtab)> asv_names <- paste0("ASV_",seq_along(asv_seqs))> asv_dna <- DNAStringSet(asv_seqs)> names(asv_dna) <- asv_names S4Vectors:::normarg_names(value, class(x), length(x)) でエラー: attempt to set too many names (1) on GroupedIRanges object of length 0> writeXStringSet(asv_dna, filepath="MAO7_dd2.fasta")
2nd try
@R% drp2 <- derepFastq('./mao7_trm_max16_cutp_2nd.fastq.gz', verbose=T)Dereplicating sequence entries in Fastq file: ./mao7_trm_max16_cutp_2nd.fastq.gzEncountered 82078 unique sequences from 98308 total sequences read.@R% err2 <- learnErrors(drp2, errorEstimationFunction=PacBioErrfun, BAND_SIZE=32, multithread=T)147899158 total bases in 98308 reads from 1 samples will be used for learning the error rates.@R% saveRDS(err2,"./mao7_err2.rds")@R% dd2 <- dada(drp2,err=err2, BAND_SIZE=32,multithread=T)Sample 1 - 98308 reads in 82078 unique sequences.@R% saveRDS(dd2,"MAO7_dd2.rds")@R% st2 <- makeSequenceTable(dd2); dim(st2)2nd try
3rd try (removed primers)
3rd try (removed primers)
2025/6/17
> drp2 <- derepFastq('./working',verbose=T)Dereplicating sequence entries in Fastq file: ./working/mao7_rmfr.fastqEncountered 47160 unique sequences from 65078 total sequences read.> err2 <- learnErrors(drp2, errorEstimationFunction=PacBioErrfun, BAND_SIZE=32, multithread=T)96635688 total bases in 65078 reads from 1 samples will be used for learning the error rates.The max qual score of 93 was not detected. Using standard error fitting. The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.The max qual score of 93 was not detected. Using standard error fitting.> saveRDS(err2, "./mao7_err2.rds")> dd2 <- dada(drp2, err=err2, BAND_SIZE=32, multithread=T)Sample 1 - 65078 reads in 47160 unique sequences.> saveRDS(dd2, "MAO7_dd2.rds")> st2 <- makeSequenceTable(dd2); dim(st2)[1] 1 104> seqtab <- readRDS("MAO7_dd2.rds")> asv_seqs <- colnames(seqtab)> asv_names <- paste0("ASV_", seq_along(asv_seqs))> asv_dna <- DNAStringSet(asv_seqs)> names(asv_dna) <- asv_namesS4Vectors:::normarg_names(value, class(x), length(x)) でエラー: attempt to set too many names (1) on GroupedIRanges object of length 0> head(asv_names)[1] "ASV_"> asv_names <- paste("ASV_", seq_along(asv_seqs))> names(asv_dna) <- asv_namesS4Vectors:::normarg_names(value, class(x), length(x)) でエラー: attempt to set too many names (1) on GroupedIRanges object of length 0> asv_dna <- DNAStringSet(asv_seqs)> names(asv_dna) <- asv_namesS4Vectors:::normarg_names(value, class(x), length(x)) でエラー: attempt to set too many names (1) on GroupedIRanges object of length 0> writeXStringSet(asv_dna, filepath="MAO7_dd2.fasta")
4th try (removed primers)
4th try (removed primers)
2025/6/18
dada2=ver. 1.36, @R
Ubuntuのupdateに伴い、Bioconductorを再インストールした。
Rのライブラリ読み込み
@R% library(dada2);library(Biostrings); library(ShortRead)
@R% library(ggplot2); library(reshape2)
@R% drp <- derepFastq("./working/",verbose=T) # 重複配列除去
Dereplicating sequence entries in Fastq file: ./working//mao7_rmfr.fastq
Encountered 47160 unique sequences from 65078 total sequences read.
@R% err <- learnErrors(drp, multithread=T) # エラーモデル学習
96635688 total bases in 65078 reads from 1 samples will be used for learning the error rates.
@R% plotErrors(err)
@R% saveRDS(err, "./mao7_err.rds")
@R% dd2 <- dada(drp, err=err, multithread=T) # ASV推定: Error correction
Sample 1 - 65078 reads in 47160 unique sequences.
@R% saveRDS(dd2, "./mao7_dd2.rds")
@R% seqtable <- makeSequenceTable(dd2);dim(seqtable)
[1] 1 104
@R% seqtable_nonchim <- removeBimeraDenovo(seqtable, method="consensus",multithread=T); dim(seqtable_nonchim)
[1] 1 17
#-- Writeout nonchimera-seq
@R% asv_names <- paste("ASV_", seq_along(seqtable_nonchim))
@R% asv_dna <- DNAStringSet(colnames(seqtable_nonchim))
@R% names(asv_dna) <- asv_names
@R% writeXStringSet(asv_dna,"MAO7_nonchim.fasta")
#-- Writeout seq (untreated chimera removal)
@R% asv_names <- paste("ASV_", seq_along(seqtable))
@R% asv_dna <- DNAStringSet(colnames(seqtable))
@R% names(asv_dna) <- asv_names
@R% writeXStringSet(asv_dna,"MAO7_onlyasv.fasta")
Page updated
Report abuse