Yersinia pestis KIM D27, using Illumina paired-end reads, with CA8.1

From wgs-assembler
Jump to: navigation, search


We will be assembling a 600bp Illumina paired-end library from Yersinia pestis KIM D27, SRP001358. This is captured in experiment SRX048908.

Fetch

mkdir READS-sra
cd READS-sra

curl -O ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR133/SRR133640/SRR133640.sra

Convert

Using fastq-dump from the SRA Toolkit, convert the NCBI .sra files into .fastq files, then generate a FRG format wrapper to load them into the assembler.

% fastq-dump \
    --split-files --split-spot --skip-technical \
    --minReadLen 64 \
    --defline-seq  @\$ac.\$si \
    --defline-qual "+" \
    ./SRR133640.sra

% fastqToCA \
    -libraryname SRR133640 \
    -insertsize 600 60 \
    -technology illumina \
    -type sanger \
    -mates SRR133640_1.fastq,SRR133640_2.fastq \
  > SRR133640.frg

Correct

Instead of using Overlap Based Trimming to trim the reads, we'll try to correct errors, trimming off ends that cannot be corrected. As with most other correction strategies, merTrim uses k-mers as evidence. We first concatenate the short reads into larger sequences (using fastq-to-fasta-merged.pl), then count k-mers, and finally dump the histogram of counts.

% cat SRR133640_1.fastq SRR133640_2.fastq | fastq-to-fasta-merged.pl > SRR133640.merged.fasta

% meryl -v -B -C -m 19 -s SRR133640.merged.fasta -o SRR133640
Have NO LIMITS!: mersPerBatch=233487221 segmentLimit=1
basesPerBatch = 233491396
Computing 1 segments using AS MUCH MEMORY AS NEEDED.
  numMersActual      = 233487221
  mersPerBatch       = 233487221
  basesPerBatch      = 233491396
  numBuckets         = 8388608 (23 bits)
  bucketPointerWidth = 28
  merDataWidth       = 15
 Allocating 417MB for mer storage (15 bits wide).
 Allocating 28MB for bucket pointer table (28 bits wide).
 Allocating 32MB for counting the size of each bucket.
 Counting mers in buckets:  186.92 Mmers -- 19.14 Mmers/second
 Creating bucket pointers.
 Releasing 32MB from counting the size of each bucket.
 Filling mers into list:    186.92 Mmers --  5.37 Mmers/second
 Writing output:            186.92 Mmers -- 19.06 Mmers/second
Segment 0 finished.

% meryl -Dh -s SRR133640 > SRR133640.histogram
Found 186922193 mers.
Found 25181862 distinct mers.
Found 19096862 unique mers.
Largest mercount is 4312; 0 mers are too big for histogram.

The histogram is showing the expected shape, with coverage peak at about 33x, and LOTS of error kmers.

K-mer histogram

Correction is done once per file, then we create a new FRG format wrapper for the corrected reads.

% merTrim -v -t 8 -m 19 -mc SRR133640 -mCillumina -F SRR133640_1.fastq -o SRR133640_corrected_1.fastq
Guessed X coverage is 34
Use minCorrect=11 minVerified=8
creating adapter mer database.
tblBits 21 seqlen 1107768
Compressed from 739092 to 5775 (13 bits)
tblBits 20 seqlen 1107768
Compressed from 739092 to 5775 (13 bits)
loading genome mer database from meryl 'SRR133640'.
 17229.6/s -      603 queued for compute;  1155296 finished;      115 queued for output)ritten;      115 queued for output)
Success!  Bye.

% merTrim -v -t 8 -m 19 -mc SRR133640 -mCillumina -F SRR133640_2.fastq -o SRR133640_corrected_2.fastq
Guessed X coverage is 34
Use minCorrect=11 minVerified=8
creating adapter mer database.
tblBits 21 seqlen 1107768
Compressed from 739092 to 5775 (13 bits)
tblBits 20 seqlen 1107768
Compressed from 739092 to 5775 (13 bits)
loading genome mer database from meryl 'SRR133640'.
 17160.9/s -      660 queued for compute;  1155239 finished;      128 queued for output)ritten;      128 queued for output)
Success!  Bye.

% fastqToCA \
    -libraryname SRR133640 \
    -insertsize 600 60 \
    -technology illumina \
    -type sanger \
    -mates SRR133640_corrected_1.fastq,SRR133640_corrected_2.fastq \
  > SRR133640.corrected.frg

cd ..

Assemble

We'll assemble both the raw uncorrected reads and the merTrim correcrted reads. The uncorrected reads are first trimmed with Overlap Based Trimming, while the corrected reads are assembled as is.

SGE is enabled, but no special settings are needed.

runCA -p ypestis -d ypestis-raw       useGrid=1 scriptOnGrid=1 doOBT=1 unitigger=bogart READS-sra/SRR133640.frg
runCA -p ypestis -d ypestis-corrected useGrid=1 scriptOnGrid=1 doOBT=0 unitigger=bogart READS-sra/SRR133640.corrected.frg

Results

Far from great assemblies, especially compared to the assembly with 454 reads.

QC statistics comparison

With corrected reads, the assembly is much cleaner (fewer small scaffolds, small contigs, degenerate contigs, singletons), has more good mate pairs incorporated, and more reads assembled into contigs. None of this helped the assembly structure though.

mergeqc.pl -wiki ypestis-*/*qc > qc.wiki
Files
ypestis-corrected/ypestis.qc ypestis-raw/ypestis.qc
Scaffolds
TotalScaffolds 145 146
TotalContigsInScaffolds 289 315
MeanContigsPerScaffold 1.99 2.16
MinContigsPerScaffold 1 1
MaxContigsPerScaffold 9 9
TotalBasesInScaffolds 4486873 4476841
MeanBasesInScaffolds 30944 30663
MinBasesInScaffolds 1655 193
MaxBasesInScaffolds 120993 120082
N25ScaffoldBases 71468 71451
N50ScaffoldBases 49329 51552
N75ScaffoldBases 31714 34140
ScaffoldAt1000000 73775 74742
ScaffoldAt2000000 52093 56764
ScaffoldAt3000000 37107 39292
ScaffoldAt4000000 17946 18052
TotalSpanOfScaffolds 4484459 4480003
MeanSpanOfScaffolds 30927 30685
MinScaffoldSpan 1655 193
MaxScaffoldSpan 120930 120053
IntraScaffoldGaps 144 169
2KbScaffolds 144 141
2KbScaffoldSpan 4482804 4475127
MeanSequenceGapLength -17 19
Top5Scaffolds=contigs,size,span,avgContig,avgGap,EUID
0 6,120993,120930,20166,-13,7180000005736 6,120082,120053,20014,-6,7180000014525
1 5,115726,115652,23145,-18,7180000005773 4,111978,111918,27994,-20,7180000014555
2 3,111924,111884,37308,-20,7180000005713 6,103338,103238,17223,-20,7180000014522
3 5,103365,103285,20673,-20,7180000005722 7,102219,102178,14603,-7,7180000014553
4 8,102377,102281,12797,-14,7180000005728 4,97479,97419,24370,-20,7180000014497
total 27,554385,554032,20533,-16 27,535096,534806,19818,-13
Contigs
TotalContigsInScaffolds 289 315
TotalBasesInScaffolds 4486873 4476841
TotalVarRecords 1827 13582
MeanContigLength 15526 14212
MinContigLength 67 75
MaxContigLength 68862 67393
N25ContigBases 42689 42616
N50ContigBases 26666 25632
N75ContigBases 13364 12762
ContigAt1000000 45079 44523
ContigAt2000000 29667 28132
ContigAt3000000 18395 17483
ContigAt4000000 8315 7392
BigContigs_greater_10000
TotalBigContigs 146 139
BigContigLength 3787278 3630127
MeanBigContigLength 25940 26116
MinBigContig 10057 10061
MaxBigContig 68862 67393
BigContigsPercentBases 84.41 81.09
SmallContigs
TotalSmallContigs 143 176
SmallContigLength 699595 846714
MeanSmallContigLength 4892 4811
MinSmallContig 67 75
MaxSmallContig 9983 9975
SmallContigsPercentBases 15.59 18.91
DegenContigs
TotalDegenContigs 702 1056
DegenContigLength 173729 211268
DegenVarRecords 317 1012
MeanDegenContigLength 247 200
MinDegenContig 65 64
MaxDegenContig 14409 14394
DegenPercentBases 3.87 4.72
Top5Contigs=reads,bases,EUID
0 31984,68862,7180000004718 29849,67393,7180000014257
1 31311,66961,7180000005563 29408,66682,7180000014391
2 29872,63712,7180000005683 30136,63845,7180000014180
3 28658,61840,7180000005527 28419,63803,7180000014310
4 28626,61745,7180000005520 27185,62123,7180000014264
total 150451,323120 144997,323846
UniqueUnitigs
TotalUUnitigs 3347 10965
MinUUnitigLength 64 64
MaxUUnitigLength 68865 35388
MeanUUnitigLength 1391 484
SDUUnitigLength 5206 2017
Surrogates
TotalSurrogates 668 1076
SurrogateInstances 3162 3381
SurrogateLength 102294 138695
SurrogateInstanceLength 469221 494084
UnPlacedSurrReadLen 6246123 6200936
PlacedSurrReadLen 1036437 790637
MinSurrogateLength 77 64
MaxSurrogateLength 5450 4850
MeanSurrogateLength 153 129
SDSurrogateLength 269 205
Mates
ReadsWithNoMate 53710(2.39%) 122954(5.73%)
ReadsWithGoodMate 2007882(89.18%) 1838638(85.73%)
ReadsWithBadShortMate 0(0.00%) 0(0.00%)
ReadsWithBadLongMate 452(0.02%) 184(0.01%)
ReadsWithSameOrientMate 194(0.01%) 118(0.01%)
ReadsWithOuttieMate 254(0.01%) 196(0.01%)
ReadsWithBothChaffMate 1562(0.07%) 1340(0.06%)
ReadsWithChaffMate 13964(0.62%) 22130(1.03%)
ReadsWithBothDegenMate 104144(4.63%) 85844(4.00%)
ReadsWithDegenMate 14394(0.64%) 14070(0.66%)
ReadsWithBothSurrMate 31746(1.41%) 39812(1.86%)
ReadsWithSurrogateMate 8986(0.40%) 6846(0.32%)
ReadsWithDiffScafMate 14104(0.63%) 12670(0.59%)
ReadsWithUnassignedMate 0(0.00%) 0(0.00%)
TotalScaffoldLinks 8 0
MeanScaffoldLinkWeight 15.62 0.00
Reads
TotalReadsInput NA NA
TotalUsableReads 2251392 2144802
AvgClearRange 98 96
ContigReads 2076770(92.24%) 1972731(91.98%)
BigContigReads 1778103(78.98%) 1623983(75.72%)
SmallContigReads 298667(13.27%) 348748(16.26%)
DegenContigReads 102260(4.54%) 93286(4.35%)
SurrogateReads 74317(3.30%) 73163(3.41%)
PlacedSurrogateReads 10633(0.47%) 8490(0.40%)
SingletonReads 8678(0.39%) 14112(0.66%)
ChaffReads 8662(0.38%) 14103(0.66%)
Coverage
ContigsOnly 45.47 42.15
Contigs_Surrogates 46.87 43.54
Contigs_Degens_Surrogates 47.28 43.49
AllReads 49.30 45.84
TotalBaseCounts
BasesCount NA NA
ClearRangeLengthFRG NA NA
ClearRangeLengthASM 221202093 205210547
SurrogateBaseLength 7282560 6991573
ContigBaseLength 204039427 188705044
DegenBaseLength 10056265 8983617
SingletonBaseLength 860278 1320950
Contig_SurrBaseLength 210285550 194905980
gcContent
Content 47.47 47.40
Unitig Consensus
NumColumnsInUnitigs 17376329 20939861
NumGapsInUnitigs 1477 4138
NumRunsOfGapsInUnitigReads 61272 286103
Contig Consensus
NumColumnsInUnitigs 4665703 4698344
NumGapsInUnitigs 5104 10103
NumRunsOfGapsInUnitigReads 244572 384187
NumColumnsInContigs 4665554 4698206
NumGapsInContigs 4953 9962
NumRunsOfGapsInContigReads 235688 374194
NumAAMismatches 2672 21377
NumVARStringsWithFlankingGaps 1141 1299