Dpuiu CA

Links

• BOG paper

Data Issues

• Problems if the coverage > read length; solution: sample the reads
• Analyze the data before assemble it:
  • If there is a reference genome map the reads to it
  • Plot kmer frequencies (use multiple kmer sizes) ; see how they compare with the estimated coverage

Sourceforge

• Wiki
• WGS Sourceforge Wiki
• WGS Sourceforge Wiki Release History

• Constants:

 cat src/AS_global.h
 ...
 #define AS_READ_MIN_LEN                   64  # should decrease for short Illumina reads  32 ???
 #define AS_OVERLAP_MIN_LEN                40  # should decrease for short Illumina reads  30 ???  

 #define AS_READ_MAX_PACKED_LEN            104  # max Illumina read length; incread to 124
 
 #define AS_READ_MAX_NORMAL_LEN_BITS       11   # should be increased to 12+ if there are many frequent kmesra nd dpn't fit in the hash

cat src/AS_CGW/AS_CGW_dataTypes.h
#define CGW_MIN_DISCRIMINATOR_UNIQUE_LENGTH 1000   # should be decraesed for short reads assemblies ; all unitigs less than that do not get scaffolded

• Download & compile CVS tip

 cvs -z3 -d:pserver:anonymous@wgs-assembler.cvs.sourceforge.net:/cvsroot/wgs-assembler co -P wgs-assembler

 cd wgs-assembler/src
 make SITE_NAME=LOCAL 

 => executables under wgs-assembler/Linux-amd64/bin/

• Download & compile release

 wget ...

 bzip2 -dc wgs-6.0-beta.tar.bz2 | tar -xf -
 cd wgs-6.0-beta
 cd kmer
 sh configure.sh
 gmake install
 cd ../src
 gmake
 cd ..

• Compile for debugging :

 make SITE_NAME=LOCAL  BUILDDEBUG=1  # in "c_make.as"

Run

• Help

runCA.pl -fields
runCA.pl -options

• Batch submission

 touch runCA.sh
   echo --- runCA $$ restarted by $USER at `date` --- >>& runCA.log
   ~/bin/runCA -d . -p prefix -s ~/bin/runCA.spec  *.frg >>& runCA.log

 at runCA.sh   

• Sanger

 runCA -d . -p prefix -s spec Sanger.frg

• doOverlapTrimming=1 (default)
• low qual Sanger reads (ex: 20 or 'D') are discarded
• Illumina reads are trimmed at 5'/3' (bp with qual < 5 (E) are deleted; read is deleted if new length< 64)
• astatLowBound 1 (default)  ; can be lowered for low/un-uniform coverage
• astatHighBound 5 (default)
• computeInsertSize=0 (default)

• 454 & Hybrid : convert sff to frg

 sffToCA -libraryname 454 -output 454 454.sff
 runCA -d . ovlOverlapper=mer unitigger=bog -p prefix *.frg
 runCA -d . -s ~/bin/runCA.454.spec -p prefix *.frg 

• gatekeeper -L : search for 454 paired end linker
• AS_GKP_sff.c  : Load_SFF function
• for mated libs set mean=3000,stdev=300

Parameters

• stopAfter:

   initialStoreBuilding
   overlapBasedTrimming | OBT     
   overlapper                         # runs overlap, overlapStore
   unitigger                          # runs frg/ovl correction
   consensusAfterUnitigger
   scaffolder
   consensusAfterScaffolder

• Example:

 ~/bin/runCA -d . -p asm -s ~/bin/runCA.454.spec stopAfter=initialStoreBuilding *.frg

Gatekeeper

Build

• Create & add

 gatekeeper    -o prefix.gkpStore -T -F prefix*frg
 gatekeeper -a -o prefix.gkpStore -T -F prefix*frg

Update

• edit frg: act:R !!!

 {DST
 act:R
 acc:10691
 mea:161662
 std:28760
 }

• append

 gatekeeper -a -o prefix.gkpStore -T -F edit.frg

• Using tab file

 Example: No restrictions on lib stddev with --edit option

 cat prefix.update
 lib     uid     19070   mean    167001
 lib     uid     19070   stddev  167001
 lib     uid     19070   isNotRandom     1
 lib     uid     19070   Orientation     O            # outie orientation for Illimina libs >1K
    
 gatekeeper --edit prefix.update prefix.gkpStore

• Update read CLR/CLV

 gatekeeper -a -r prefix.clr -o prefix.gkpStore   # not in the wgs release
 gatekeeper -a -v prefix.clv -o prefix.gkpStore

Delete

• DST messages can not be deleted
• A FRG can be deleted only if there is no related LKG message; an abbreviated MSG
• act:D !!!

 # frg1
 {LKG
 act:D
 typ:M
 fg1:94376
 fg2:94654
 }
 {FRG
 act:D
 acc:94376
 }

 # frg2 ???
 {LKG
 act:D
 fg:94376
 fg:94654
 }
 {FRG
 act:D
 acc:94376
 }

 cat delete.ids | perl -ane 'print "{FRG\nact:D\nacc:$F[0]\n}\n";'> delete.frg
 cat delete.ids | ~/bin/deleteFrg.pl -forward 1 -reverse 2 > delete.frg            # forward reads have suffix 1; forward reads have suffix 1

 gatekeeper -a -o prefix.gkpStore -T -F delete.frg

Dump

• INFO

 gatekeeper -dumpinfo prefix.gkpStore

• LAST frg in store

 gatekeeper -dumpinfo -lastfragiid prefix.gkpStore

• FASTA sequences

 gatekeeper -dumpfastaseq prefix.gkpStore > prefix.seq
 gatekeeper -dumpfastaqlt prefix.gkpStore > prefix.qual

• FASTA reads for newbler

 gatekeeper -dumpnewbler newbler prefix.gkpStore => prefix.fna,prefix.fna.qual

 grep "^>" prefix.fna | head -2
 >SRR001355.3635.1a template=2020+2021 dir=F library=SRX000348 trim=20-117
 >SRR001355.3635.1b template=2020+2021 dir=R library=SRX000348 trim=1-130

• FASTQ reads

 gatekeeper -dumpfastq prefix. prefix.gkpStore
 =>
 prefix.1.fastq
 prefix.2.fastq
 prefix.paired.fastq    # contains both asm.1.fastq & asm.2.fastq reads
 prefix.unmated.fastq

 cat prefix.1.fastq | perl ~/bin/clrFastq.pl > prefix.1.clr.fastq 
 cat prefix.2.fastq | perl ~/bin/clrFastq.pl > prefix.2.clr.fastq 
 cat prefix.paired.fastq | perl ~/bin/clrFastq.pl > prefix.paired.clr.fastq 
 cat prefix.unmated.fastq | perl ~/bin/clrFastq.pl > prefix.unmated.clr.fastq 

• FRG file

 gatekeeper -dumpfrg prefix.gkpStore > prefix.frg
 gatekeeper -dumpfrg prefix.gkpStore -format2 > prefix.frg2 

• READS in TAB format

 gatekeeper -dumpfragments -tabular prefix.gkpStore > prefix.gkpStore.info
 gatekeeper -dumpfragments -tabular prefix.gkpStore | grep -v ^UID | awk '{print $1,$2}' > prefix.gkpStore.uid2iid

• LIBS in TAB format

 gatekeeper -dumplibraries -tabular prefix.gkpStore | grep -v ^UID | awk '{print $1,$4,$5}' | sed 's/\.000//g'> prefix.dst  # original vals
 asm2mdi.pl < prefix.asm > prefix.mdi                                                                                       # final vals

• READ CLR's (default=LATEST)

 gatekeeper -dumpfragments -tabular prefix.gkpStore/ | perl -ane 'print join "\t", @F[1,0,-2,-1,6]; print "\n";' | head | pretty -o
 IID    UID                                  clrBeginLATEST  clrEndLATEST  isDeleted  
 1      HWI-EAS385_0001:1:1:1569:15197#0/1   0               65            0          
 2      HWI-EAS385_0001:1:1:2095:12498#0/1   0               67            0          
 3      HWI-EAS385_0001:1:1:5452:4010#0/1    0               78            0          
 4      HWI-EAS385_0001:1:1:12400:16773#0/1  0               105           0          
 5      HWI-EAS385_0001:1:1:12884:18158#0/1  0               80            0          
 ...
 
 gatekeeper -dumpfragments -tabular prefix.gkpStore/  | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > prefix.clr
 gatekeeper -dumpfragments -tabular -clear [CLR|OBTINITIAL|OBTMERGE|OBTCHIMERA|ECR_0|ECR_1|ECR_2|LATEST] prefix.gkpStore/  | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > prefix.?.clr

 gatekeeper -dumpfragments -tabular -clear CLR asm.gkpStore        | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > asm.CLR.clr
 gatekeeper -dumpfragments -tabular -clear OBTINITIAL asm.gkpStore | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > asm.OBTINITIAL.clr &
 gatekeeper -dumpfragments -tabular -clear OBTCHIMERA asm.gkpStore | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > asm.OBTCHIMERA.clr &
 ...
 gatekeeper -dumpfragments -tabular -clear LATEST asm.gkpStore     | perl -ane 'print join "\t", @F[0,-2,-1]; print "\n";' > asm.LATEST.clr &

• READ CLR's summary

 gatekeeper -dumpfragments asm.gkpStore | grep ^fragmentClear | perl -ane '  print $1," ",$3-$2,"\n" if(/(\w+),(\d+),(\d+)$/);' | ~/bin/sum2.pl 

• MATED

 frg22mates.pl < prefix.frg > prerfix.mates

• Deleted reads

 gdt asm.gkpStore | p 'print $F[0],"\n" if($F[6]);'

• All scaffold reads in FRG2 format

 cat asm.posmap.frgscf | grep scfid > scfid.posmap.frgscf
 gatekeeper -uid scfid.posmap.frgscf -dumpfrg -donotfixmates -format2 asm.gkpStore > scfid.frg

Id mapping

• Frg FASTQ packed format id's are no recorded in the gatekeeper; assgned a new UID ; to map them back ...

 ~/bin/fastq2wgsUIDs.pl s_*txt > s.uid2iid2id
 gatekeeper -dumpfragments -tabular asm.gkpStore | awk '{print $1,$2,$10}' > s.uid2iid

 paste s.uid2iid2id s.uid2iid | p 'print "ERROR $_" if($F[3] ne $F[-1]);' # check the lengths

Meryl

 # Get mer frequencies
 meryl -C -B -m 22 -s prefix.seq -o prefix.22mers 

 # Dump mer sequences
 meryl -Dt -s prefix.22mers > prefix.22mers.seq

 # Dump histogram
 meryl -Dh -s prefix.22mers | sort -nk2 -r  >  prefix.22mers.histo

 # Get missing mers
 meryl -C -B -m 22 -s ref.fasta -o ref.22mers
 meryl -C -B -m 22 -s qry.fasta -o qry.22mers
 meryl -M sub -s ref.22mers -s qry.22mers -o diff.22mers
 meryl -Dt -s diff.22mers  > diff.22mers.seq

 # Combine mers

• jellyfish

 # Get mer frequencies
 jellyfish count -s 1000000000 -m 18 --both-strands -t 16 -o prefix.jf prefix*.fastq [prefix*.seq]
 jellyfish merge prefix.jf_? -o prefix.jf

 # Dump histogram
 jellyfish histo prefix.jf 

• countKmers.pl (other)

 cat lib.seq | ~/bin/countKmers.pl -k 22 | head -20 |          grep -v "^>" > 22mers.fwd.grep
 cat lib.seq | ~/bin/countKmers.pl -k 22 | head -20 | revseq | grep -v "^>" > 22mers.rev.grep
 cat 22mers.fwd.grep 22mers.rev.grep > 22mers.grep
 cat lib.seq | ~/bin/fgrep.pl -f 22mer.grep | count.pl

---

meryl -C -B -m 12 -s prefix.seq -o prefix.12mers

meryl -Dh -s prefix.12mers | sort -nk2 -r | more

OBT

• runCA spec

 doOverlapTrimming = 1

• command

 overlap -G ...   # -G to compute partial overlaps

• Dump

 overlapStore -d ./0-overlaptrim/asm.obtStore

• Number of jobs: should be < fileListMax (1024*10 by default) otherwise overlapStore fails

 ~/bin/getOvlJobsCount.pl ovlHashBlockSize ovlRefBlockSize frgCount

 #or
 cat 0-overlaptrim-overlap/ovlopts.pl | grep -c '^"h'

Overlap

• Jobs

 merOverlapperSeedBatchSize     100000
 merOverlapperExtendBatchSize   75000
 
 overmerry.sh :       #frgs/merOverlapperSeedBatchSize   jobs -> seeds/ dir
 olap-from-seeds.sh : #frgs/merOverlapperExtendBatchSize jobs -> olaps/ dir

• Num ber of jobs

 cat 1-overlapper/ovlopts.pl | grep -c '^"h'

• Restart in runCA

 rm 1-overlapper/o*

OverlapStore

• Flags:

 -t -M            : don't improve performance much
 -I ignore_file   : not implemented

• Dump

 overlapStore -d prefix.ovlStore > prefix.ovl.tab

 overlapStore -d prefix.ovlStore | awk '{print $1}' | uniq -c | awk '{print $2,$1}' > prefix.ovlStore.count.tmp
 cat prefix.ovlStore.count.tmp | perl -ane 'foreach ($p+1..$F[0]-1) { print "$_ 0\n" } $p=$F[0]; print $_'  > prefix.ovlStore.count
 rm prefix.ovlStore.count.tmp

• Get stats; median number of ovls (5'/3') for each lib

 overlapStats -G prefix.gkpStore -O prefix.ovlStore/ -o prefix
 cat prefix.repeatmodel.lib.00*stats | egrep '^Lib|nSamples|median' | p 'chomp; if(/Lib/) {print "\n$_"} else { print "  $F[1]"}' | pretty -o

• Delete overlaps from *.obt.gz files

 zcat *.obt.gz  | convertOverlap -a -ovl |  ~/bin/difference1or2.pl -f delete.acc | convertOverlap -b -ovl | gzip > new.obt.gz

• Delete overlaps from ovlStore : asm=>new

 overlapStore -d asm.ovlStore | perl -ane 'print $_ if($F[0]<$F[1]);' | ~/bin/difference1or2.pl -f delete.iid |  convertOverlap -b -ovl | gzip  >! new.ovb.gz
 overlapStore -c new.ovlStore -g asm.gkpStore/ new.ovb.gz

 #not necessary to gzip the overlap files
 overlapStore -d asm.ovlStore | perl -ane 'print $_ if($F[0]<$F[1]);' | ~/bin/difference1or2.pl -f delete.iid |  convertOverlap -b -ovl >! new.ovb
 overlapStore -c new.ovlStore -g asm.gkpStore/ new.ovb

• Edit(invert) overlaps from ovlStore : asm=>new

 ... > invert.iid     # get the iids of the reads to change          
 overlapStore -d asm.ovlStore  | p 'print $_ if($F[0]<=$F[1]);' | ~/bin/invertOvl.pl -f invert.iid | convertOverlap -b -ovl | gzip > new.ovb.gz
 overlapStore -c new.ovlStore -g asm.gkpStore new.ovb.gz

• Count how many overlaps have errors

 overlapStore -d asm.ovlStore | p 'print if($F[5]>0);' | wc -l

• Get 5'/3' overlap stats

 overlapStore -d asm.ovlStore | perl /nfshomes/dpuiu/Archives/JCVI/bin/countOverlaps53.pl | getSummary.pl -i 1 -t 5 
 overlapStore -d asm.ovlStore | perl /nfshomes/dpuiu/Archives/JCVI/bin/countOverlaps53.pl | getSummary.pl -i 2 -t 3

• Build faster ...

 Example (Bee):
 e 'foreach ("0000000001".."0000000058") { print "overlapStore -c $_.ovlStore -g asm.gkpStore/ -i 0 -M 8192 1-overlapper/$_/*gz\n" }' | scheduler.pl  # 13min/job  
 e 'foreach ("0000000002".."0000000058") { print "overlapStore -m 0000000001.ovlStore $_.ovlStore\n" }'                                               # 3 min/job 
 total=2*13+57*3=3.5 hrs (instead of 13.85 reported in runCA.log)

Frgcorr

• Fails on deleted reads with overlaps ; must delete the overlaps as well

Unitiger

• Unitig stats:

 cat 4-unitigger/*.cga.0 | perl ~/bin/utgLog2utgLen.pl | getSummary.pl -i 0 -j 1

• Library insert estimates:

 cat 4-unitigger/unitigger.err | grepi -A 9 LIB | pretty   # in this case 9 is the number of libs

 should compare with the original values
 gatekeeper -dumplibraries -tabular *.gkpStore | more

• Bog: don't let small inserts break unitigs:

 edit AS_BOG_MateChecker.cc

 void
 MateLocation::buildHappinessGraphs(Unitig *utg, DistanceCompute *globalStats) {

 for (uint32 mleidx=0; mleidx<_table.size(); mleidx++) {
   MateLocationEntry &loc = _table[mleidx];

   uint32 lib =  _fi->libraryIID(loc.mleFrgID1);

   //Aelksey Zimin also do not use short Illumina mate pairs
   //  if (globalStats[lib].samples < 10 )
   if (globalStats[lib].samples < 10 || globalStats[lib].stddev < 100 || globalStats[lib].mean < 500)

Tigstore

• list unitigs

 tigStore -g asm.gkpStore -t asm.tigStore 2 -D unitiglist     # ids  
 maID    isPresent  isDeleted  ptID  svID  fileOffset  covStat  
 0       1          0          1     2     0           3334     
 1       1          0          1     2     521908      -1       
 2       1          0          1     2     522764      543      
 ...

 tigStore -g asm.gkpStore -t asm.tigStore 2 -D properties     # exit: stdout    # properties
 unitig_coverage_stat        0 660
 unitig_microhet_prob        0 0.000000
 unitig_status               0 X
 unitig_unique_rept          0 X
 ...

• number of tigs

tigStore -g asm.gkpStore -t asm.tigStore 2 -D unitiglist | tail -1 | awk '{print $1}'

• list frags/consensus/layout in a unitig

 tigStore -g asm.gkpStore -t asm.tigStore 2 -u 0 -d frags 
 tigStore -g asm.gkpStore -t asm.tigStore 2 -u 0 -d consensus
 tigStore -g asm.gkpStore -t asm.tigStore 2 -u 0 -d layout
 tigStore -g asm.gkpStore -t asm.tigStore 2 -u 0  -d properties  # exit: stderr & stdout

• dump all unitig consensus

 tigStore -g asm.gkpStore -t asm.tigStore 2 -U -d consensus > utg.fasta

• get single read unitigs

 tigStore -g asm.gkpStore -t asm.tigStore 43 -D properties | grep ^unitig_num_frags | grep ' 1$' 

• get read/utg or ctg stats

 tigStore -g asm.gkpStore -t asm.tigStore 43 -D properties |  grep ^unitig_num_frags | getSummary.pl -z 1 -i 2 -t unitig_num_frags 
 tigStore -g asm.gkpStore -t asm.tigStore 43 -D properties |  grep ^contig_num_frags | getSummary.pl -z 1 -i 2 -t contig_num_frags 

• delete a unitig : get & update layout

 tigStore -g asm.gkpStore -t asm.tigStore 2 -u id -d layout | tee utg.layout | grep ^unitig | awk '{print $2}' | ~/bin/deleteUnitigIdLayout.pl > utg.layout.delete
 # or
 tigStore -g asm.gkpStore -t asm.tigStore 2 -u id -d layout | tee utg.layout |  ~/bin/deleteUnitigLayout.pl > utg.layout.delete
 tigStore -g asm.gkpStore -t asm.tigStore 2 -R utg.layout.delete
 tigStore -g asm.gkpStore -t asm.tigStore 2 -D unitiglist | p 'print if($F[2]);'

• delete single frg unitigs

 tigStore -g asm.gkpStore -t asm.tigStore 2 -U -d layout | egrep '^unitig|^data.num_frags' | p 'chomp; print $_," "; print "\n" if(/^data/)' | p 'print $F[1],"\n" if($F[3]==1);' | ~/bin/deleteUnitigIdLayout.pl > utg.layout.delete 
 tigStore -g asm.gkpStore -t asm.tigStore 2 -R utg.layout.delete

• assign unitig -1 to all unitigs (for addition to a new store)

 tigStore -g asm.gkpStore -t asm.tigStore 2 -u id -d layout | tee utg.layout | p 'if (/^unitig/) { print "unitig -1\n"} elsif(/^UTG/) { $F[7]=-1; print join " ",@F; print "\n" } else { print $_};' > utg.layout.add
 tigStore -g asm.gkpStore -t new.tigStore 2 -R utg.layout.add

 tigStore -g asm.gkpStore -t asm.tigStore 1 -u id -d layout | tee utg.layout | p 'if (/^unitig/) { print "unitig -1\n"} else { print $_};' > utg.layout.add
 tigStore -g asm.gkpStore -t new.tigStore 1 -R utg.layout.add

• fix unitig consensus

 touch utg.layout
 grep FAILED 5-consensus/asm*err | cut -f3 -d ' ' | sort -u -n | p 'system "tigStore -g asm.gkpStore -t asm.tigStore 2 -u $F[0] -d layout >> utg.layout\n";'  
 cat  utg.layout | updateUTGlayout.pl  >  utg.layout.update
 tigStore -g asm.gkpStore -t asm.tigStore 2 -R utg.layout.update

 grep FAILED 5-consensus/asm*err | sed 's/\.err:/ /' | cut -f1 | sort -u | p 'system "touch $F[0].success\n";' 
 touch 5-consensus/consensus.success

• Filter 2+ fragment unitigs

CGW

• Long run times usually mean problems (conflicts) with the data
• Drbug

 gdb /fs/szdevel/dpuiu/SourceForge/wgs/Linux-amd64/bin//cgw
 $ run -j 1 -k 5 -r 5 -s 2 -z -P 2 -B 75000 -m 100 -g ./asm.gkpStore -t ./asm.tigStore -o asm
 bt

• Differences between the LAST run of CGW and PREVIOUS ones

 # options
 PREV       LAST
            -R 7             # restart from checkpoint file number 'ckp'
            -N ckp01-ABS     # restart from checkpoint location 'ckp' (see the timing file)
 -s 0       -s 2             # stone throwing level
 -S 0                        # S==0 : do not resolve surrogates
 -G                          # Don't generate output (cgw or cam) 

 #files
 difference.pl 7-4-CGW.ls 7-0-CGW.ls
 asm.asm.cam
 asm.dregs.cam
 asm.partitionInfo  #  "equivalent of" 4-unitigger/asm.partitioningInfo 
 asm.partitioning   #  "equivalent of" 4-unitigger/asm.partitioning

• running only one instance of ECR

 rm -fr 7-[234]*
 /fs/szdevel/dpuiu/SourceForge/wgs-assembler.030210/Linux-amd64/bin/runCA -d . -p asm -s ./runCA.spec doExtendClearRanges=1 *.frg >> & runCA.log

• canceling ECR run (if doExtendClearRanges set to 1 or 2 and the run takes too long)

 grep ckp ./7-0-CGW/cgw.out 
 ./7-0-CGW/asm.ckp.3 (logical ckp01-ABS) after building scaffolds at ...
 ./7-0-CGW/asm.ckp.4 (logical ckp03-ACD) after scaffold cleaning at ...
 ./7-0-CGW/asm.ckp.5 (logical ckp05-1SM) after 1st scaffold merge at ...
 ./7-0-CGW/asm.ckp.6 (logical ckp06-AS) after stone throwing at ...
 ./7-0-CGW/asm.ckp.7 (logical ckp08-2SM) after 2nd scaffold merge at ...
 ./7-0-CGW/asm.ckp.8 (logical ckp09-FR) after final rocks at ...

 cgw -R 7 -N ckp09-FR -j 1 -k 5 -r 5 -s 2 -z -P 2 -B 75000 -m 100 -g ./asm.gkpStore -t ./asm.tigStore -o ./7-0-CGW/asm

 grep ckp ./7-0-CGW/cgw.out 
 ...
 ./7-0-CGW/asm.ckp.9 (logical ckp10-PS) after partial stones at ...
 ./7-0-CGW/asm.ckp.10 (logical ckp11-FCS) after final contained stones at ...
 ./7-0-CGW/asm.ckp.11 (logical ckp13-RS) after resolve surrogates at ...
 ./7-0-CGW/asm.ckp.12 (logical ckp14-FIN) after output at ...

• running ECR (if initially doExtendClearRanges=0) ; make sure cgwPurgeCheckpoints=0

 runCA -d . -p asm -s runCA.spec stopAfter=scaffolder cgwPurgeCheckpoints=0 doExtendClearRanges=0 *.frg >> & runCA.log
 
 rm 7-CGW
 cd 7-0-CGW
 mkdir OLD
 mv asm.ckp.9 asm.ckp.1? OLD asm.*.cam asm.partition* asm.timing OLD
 cd ..
 
 runCA -d . -p asm -s runCA.spec stopAfter=scaffolder cgwPurgeCheckpoints=0 doExtendClearRanges=1 *.frg >> & runCA.log

• DistUpadte file

 7-0-CGW/stat/scaffold_final.distupdate.dst

• Stats

 cat 7-0-CGW/rezlog/stone.i02.log | perl ~/bin/cgwLog2scfLen.pl | getSummary.pl -t scfLen

ECR

• doExtendClearRanges = 2 => cgw run 3 times, ECR 2 times

Consensus

• Dump unitig seq

 tigStore -g asm.gkpStore -t asm.tigStore/ 20 -D unitigs > unitigs.txt
 cat unitigs.txt | grep -v ^maID | perl -ane 'system "tigStore -g asm.gkpStore -t asm.tigStore 20 -d consensus -u  $F[0]\n";' >& unitigs.fasta.tmp
 cat unitigs.fasta.tmp  | grep -v null | sed 's/cns=//' | sed 's/-//g' | nl0 | tab2fasta.pl > unitigs.fasta
 rm unitigs.fasta.tmp

• Dump contig seq

 tigStore -g asm.gkpStore -t asm.tigStore/ 20 -D contigs > contigs.txt
 cat contigs.txt | grep -v ^maID | perl -ane 'system "tigStore -g asm.gkpStore -t asm.tigStore 20 -d consensus -c  $F[0]\n";'>& contigs.fasta.tmp
 cat contigs.fasta.tmp | grep -v null | sed 's/cns=//' | sed 's/-//g' | nl0 | tab2fasta.pl > contigs.fasta
 rm contigs.fasta.tmp

• Recompute unitig consensus

 utgcns -g ./asm.gkpStore -t ./asm.tigStore 2 011 -u 122533 -v -V

Log file & runtime (wgs-5.2)

• Filter log

 ~/bin/getCAruntimes.pl -filter runCA.log > runCA.filter.log

• Display log

 cat runCA.log | ~/bin/getCAruntimes.pl -hour -total -plot  # => runCA.runTimes ; runCA.runTimes.gp ; runCA.runTimes.png
 display runCA.runTimes.png

• runCA.log runCA.times bos taurus BCM.CLONEEND 16,875 Sanger reads
• runCA.log runCA.times wolbachia symbioint of culex pipiens quinquefasciatus 36,767 Sanger reads
• runCA.runTimes.png Bumblebee assembly

Output processing

Posmap

• Get posMap files

buildPosMap -o prefix < prefix.asm

• Get number of reads/ctg or scf (from the pos files)

 cat prefix.posmap.frgctg | count.pl -i 1 > prefix.posmap.ctgfrgcnt
 cat prefix.posmap.frgscf | count.pl -i 1 > prefix.posmap.scffrgcnt
 cat prefix.posmap.ctgscf | count.pl -i 1 > prefix.posmap.scfctgcnt

• Compute read/mate/ctg cvg (from the pos files)

 posmap2cvg.pl                                  < prefix.posmap.frgscf > prefix.posmap.scffrgcvg
 OLD/posmap2cvg.4.pl -mates prefix.posmap.mates < prefix.posmap.frgscf > prefix.posmap.scfmatecvg
 posmap2cvg.pl                                  < prefix.posmap.ctgscf > prefix.posmap.scfctgcvg

• Lib assembly stats *

 cd 9-*
 gatekeeper -dumplibraries -tabular ../genome.gkpStore |                    pretty > genome.gkpStore.libinfo1
 gatekeeper -dumpinfo               ../genome.gkpStore | grep -A 100 ^Lib | pretty > genome.gkpStore.libinfo2

 gdt genome.gkpStore | p 'print "$F[0] $F[5]\n" unless($F[6]);' > genome.gkpStore.frg2lib
 paste genome.gkpStore.frg2lib genome.posmap.frags | p ' print "$F[1] $F[5]\n";' | count.pl >  genome.gkpStore.frgcount
 gdt genome.gkpStore | p 'print "$F[4] del\n" if($F[6] eq "1");' | count.pl >> genome.gkpStore.frgcount
 cat genome.gkpStore.frgcount | count2col.pl -i 0 -j 1 -k 2 | awk '{print $1,$4,$2,$3}' | pretty > genome.gkpStore.libcount 

 more genome.gkpStore.lib*  

• Posmap2gff

 #getting frglib info
 gdt ../genome.gkpStore | p 'print "$F[0] $F[5]\n" unless($F[6]);' > genome.gkpStore.frg2lib

 # one of ...
 cat genome.posmap.frgscf | ~/bin/posmap2gff.pl -sou frg -fea frg -p > genome.posmap.frgscf.gff
 cat genome.posmap.frgscf | ~/bin/posmap2gff.pl -sf genome.gkpStore.frg2lib -fea frg -p | sed 's/[ab]//' > genome.posmap.frgscf.gff
 cat genome.posmap.frgscf | ~/bin/posmap2gff.pl -sf genome.gkpStore.frg2lib -fea frg -p >! genome.posmap.frgscf.gff

 cat genome.posmap.utgscf | ~/bin/posmap2gff.pl -sou utg -p >! genome.posmap.utgscf.gff
 cat genome.posmap.ctgscf | ~/bin/posmap2gff.pl -sou ctg -p >! genome.posmap.ctgscf.gff

 cat genome.posmap.scflen | perl -ane 'print "$F[0]\t$F[0]\t0\t$F[1]\tf\n";' | posmap2gff.pl -sou scf >>  genome.posmap.scflen.gff

 setenv SEQID 1374988
 cat genome*gff | p 'print $_ if($F[8]==$ENV{SEQID});' >! $SEQID.gff

 ~/szdevel/Apollo/apollo -i gff -f $SEQID.gff

• Get 0cvg regions

 posmap2cvg.pl -f prefix.posmap.scflen -Max 0 < prefix.posmap.frgscf  # scaff gaps
 posmap2cvg.pl -f prefix.posmap.ctglen -Max 0 < prefix.posmap.frgctg  # surrogates

• Get links

 ~/bin/scf2scflnk.sh prefix                      
 ~/bin/ctg2ctglnk.sh prefix                      
 ...

• Add ctg/scf into to the mates fils

 join12.pl prefix.posmap.mates prefix.posmap.frgctg > prefix.posmap.mates2.tmp
 join12.pl prefix.posmap.mates2.tmp prefix.posmap.frgscf > prefix.posmap.mates2
 rm prefix.posmap.mates2.tmp

• Display graphs

cat prefix.posmap.scf2scflnk | ~/bin/tab2dot.pl | dot -Tpng -o prefix.posmap.scf2scflnk.png
cat prefix.posmap.ctg2ctglnk | ~/bin/tab2dot.pl | dot -Tpng -o prefix.posmap.ctg2ctglnk.png

• Qc stats based on posmap files

 posmap2qc.pl -ctgscf prefix.posmap.ctgscf > prefix.qc

• Break scf/ctg pipeline

 ~/bin/breakPosmapKeep.amos

• Get scaff file

 ~/bin/posmap2scaff.pl prefix.posmap.ctgscf

• Gaps longer than a certain size

 cat asm.posmap.ctgscf | p 'print $p,$_,"\n" if($F[1] eq $P[1] and $F[2]-$P[3]>1000); $p=$_;@P=@F;'

Mate redundancy

• Sequence based (first 32bp in fwd/rev Illumina reads; max 1 mismatch)

 paste s_2_?_8kb*.txt | perl -ane 'chomp; if($.%4==1) { s/\@(\S+)[12]//; printf("%-45s",$1); } elsif($.%4==2) { print substr($F[0],0,32),"\t",substr($F[1],0,32),"\n" } ' | sort -k2,3 | ./getDuplicates.pl > s_2_8k.tab
 cat s_2_8k.tab | perl -ane 'print $F[0],"1\n";' >  s_2_1_8k.redundant
 cat s_2_8k.tab | perl -ane 'print $F[0],"2\n";' >  s_2_2_8k.redundant

• overlap(Store) based

 #Example
 getMateRedundancy.amos -D BEGIN=1 -D END=219  -D GKPSTORE=../asm.gkpStore -D OVLSTORE=../asm.ovlStore s_2_3kb
 #=> asm.redundant_mate_pairs , asm.redundant_mate_clusters

• Alignment based

 cat prefix*.soap2 | ~/bin/AMOS/soap2sameRef.pl -pos | sort -nk2 -nk3 | p 'print $_ if(@P and $F[1]-$P[1]<=1 and abs($F[2]-$P[2])<=1); @P=@F;' > prefix.redundant

• Posmap file based

 cat asm.posmap.frgscf |  ~/bin/posmap2ovl.pl -p 75 | grep -v ^ref | ~/bin/getMateRedundancy12.pl > asm.redundant_mate_pairs 
 cat asm.redundant_mate_pairs |  perl ~/bin/cluster.pl | sort > asm.redundant_mate_clusters

Fasta

• Align deg to scaff

 nucmer -mum -l 40 -c 250 scf.fasta deg.fasta  -p scf-deg
 delta-filter -q scf-deg.delta > scf-deg.filter-q.delta
 ~/bin/delta2pairCvg+.pl < scf-deg.filter-q.delta > scf-deg.filter-q.cvg.pairs
 ~/bin/max2.pl -i 4 -j 6 < scf-deg.filter-q.cvg.pairs > scf-deg.filter-q.maxCvg.pairs
 cat scf-deg.filter-q.maxCvg.pairs | p 'print $_ if($F[5]-$F[6]<2000);' > scf-deg.filter-q.variants.pairs
 awk '{print $5,$6}' scf-deg.filter-q.variants.pairs > scf-deg.filter-q.variants.ids

• Indexing

 # contigs
 ~/bin/formatdb.pl < prefix.ctg.fasta > prefix.ctg.fasta.offset
 ~/bin/fastacmd.pl -o prefix.ctg.fasta.offset -f filter.ctg.ids < prefix.ctg.fasta > filter.ctg.fasta

 # same for scaff, deg ...

• Contaminant search : nucmer against

 /nfshomes/dpuiu/db/Ecoli.all
 /nfshomes/dpuiu/db/UniVec_Core
 /nfshomes/dpuiu/db/OtherVec

Asm

• Get message counts

cat asm.asm | grep "^{" | egrep -v 'MPS|UPS|VAR|CTP' | uniq -c | more
     1 {MDI
273375 {AFG  # reads
 36258 {AMP  # links
  2185 {UTG  # MPS+           messages
  4876 {ULK
  1443 {CCO  # VAR*,MPS+,UPS+ messages
  3277 {CLK
    15 {SCF
    40 {SLK  # CTP+           messages

• Get output seq

 asmOutputFasta -p prefix < prefix.asm

• Get output seq clr

 ~/bin/asm2clrs.pl < prefix.asm > < prefix.asm.clr

• Get new library insert sizes:

 asm2mdi.pl < prefix.asm > prefix.mdi

• Get surrogate ids

 #most of them have negative aStat's
 $ cat *.asm | grep -B 8 ^sta:S | grep ^acc | perl -ane '/(\d+)/; print "utg".$1,"\n";'

Qc

 caqc.pl -euid  prefix.asm

Scaffolds

• Get scaff files

 cavalidate -e 10 prefix
 bank2scaff prefix.bnk/ > prefix.scaff

 ~/bin/ca2scaff -i prefix.asm -o prefix.scaff

 cat prefix.scaff | grep "^>" | sed 's/>//' | awk '{print $1,$2}' > prefix.scfctg
 cat prefix.scaff | grep "^>" | sed 's/>//' | awk '{print $1,$3}' > prefix.scflen  
 cat prefix.scaff | grep "^>" | sed 's/>//' | awk '{print $1,$4}' > prefix.scfspan  

 tac prefix.scaff | p '$count=0 if(/^>/); $count++; print $_ if($count>2);' | awk '{print $4}'   # sequence gaps

AGP

 asmToAGP.pl < prefix.asm > prefix.agp

        /--- buildPosMap ------> .posmap
       /
  .asm /---- asmToAGP ------------------------------------------------> .AGP 
      \                                                              /   
       \---- cavalidate, bank2scaff?--> .scaff --- scaff2agp.pl --- 

• Other scripts:

 bank2scaff : incorrect gap sizes

 valiadteAgp.pl < prefix.agp

 infoseq 9-terminator/prefix.ctg.fasta | sed 's/ctg//' >! prefix.ctg.infoseq
 scaff2agp.pl -i prefix.ctg.infoseq < prefix.scaff > prefix.agp

Moving an assembly

• Minimum file set

 scp runCA.sh
 scp runCA.spec

 scp asm.gkpStore/*

 scp 3-overlapcorrection/asm.erates.updated 

 scp 4-unitigger/asm.iidmap
 scp 4-unitigger/asm.partitioning
 scp 4-unitigger/asm.partitioningInfo
 scp 4-unitigger/unitigger.success

 scp 5-consensus/consensus.success

 scp 7-0-CGW/asm.timing
 scp 7-0-CGW/asm.ckp.11  # last one
 scp 7-0-CGW/asm.partitioning
 scp 7-0-CGW/asm.partitionInfo
 scp 7-0-CGW/cgw.out.00.ckp.11
 scp 7-0-CGW/cgw.out

 scp asm.tigStore/*

Add reads to an assembly

 ln -s ../wgs-prev/prefix.frg .
 ln -s ../Data/new.frg

 cp ../wgs-prev/runCA.spec ../wgs-prev/runCA.sh   .

 ln -s ../wgs-prev/0-mercounts

 cp -R ../wgs-prev/prefix.gkpStore
 gatekeeper -a -o prefix.gkpStore new.frg

 ll -d ../wgs-prev/1-overlapper/000000000*                  # if multiple jobs ...
 mkdir -p 1-overlapper/0000000001                           # mkdir -p 1-overlapper/0000000002 ...  ???
 
 cat ../wgs-prev/1-overlapper/*pl | grep '^"h' | tail -1    # get last qry/ref block  => "h0056000001r0050000000"

 cd 1-overlapper/0000000001
 ln -s ../../wgs-prev/1-overlapper/0000000001/*ovb.gz .
 rm h0056000001*.ovb.gz *r0050000000.ovb.gz                 # this should change from case to case 
 cd ../..

Removing bad mates

 cat ../wgs-prev/9-terminator/asm.posmap.frags | grep bad > asm.posmap.frags.bad

 gatekeeper -dumpfragments -tabular -uid asm.posmap.frags.bad ./wgs-prev/asm.gkpStore |  grep -v ^UID | ~/bin/deleteLkg.pl -j 2 > asm.posmap.frags.bad.delete.lkg

 ln -s ../wgs-prev/*.frg 
 gatekeeper -o asm.gkpStore *.frg asm.posmap.frags.bad.delete.lkg 

 ln -s ../wgs-prev/0-mercounts/
 ln -s ../wgs-prev/asm.ovlStore/
 ln -s ../wgs-prev/3-overlapcorrection/
 ln -s ../wgs-prev/4-unitigger/
 cp -R ../wgs-prev/asm.tigStore/ .
 ln -s ../wgs-prev/5-consensus

 runCA ...

Test datasets

Porphyromonas_gingivalis_W83 : Sanger + 454

• TA  : ftp://ftp.ncbi.nih.gov/pub/TraceDB/porphyromonas_gingivalis_w83

 wget ftp://ftp.ncbi.nih.gov/pub/TraceDB/porphyromonas_gingivalis_w83/anc.porphyromonas_gingivalis_w83.001.gz
 ...
 wget ftp://ftp.ncbi.nih.gov/pub/TraceDB/porphyromonas_gingivalis_w83/xml.porphyromonas_gingivalis_w83.001.gz

• SRA : http://www.ncbi.nlm.nih.gov/sites/entrez?db=genomeprj&cmd=Link&LinkName=genomeprj_sra&from_uid=29961
• CBCB

 /fs/sztmpscratch/dpuiu/porphyromonas_gingivalis_w83/

• Complete genome at NCBI : NC_002950.2 ( 2,343,476bp 48.29%GC)

 bp_fetch.pl net::genbank:NC_002950.2 > porphyromonas_gingivalis_w83.1con

Input formatting

• https://sourceforge.net/apps/mediawiki/wgs-assembler/index.php?title=Formatting_Inputs

Read counts

 zcat anc.*.gz | ~/bin/traceanc2anc.pl SEQ_LIB_ID INSERT_SIZE  INSERT_STDEV SVECTOR_CODE | count.pl 

TA formatting

• frg FORMAT2

• TIs are frg uids

• All libs

 tracedb-to-frg.pl -xml xml.porphyromonas_gingivalis_w83.001.gz
 tracedb-to-frg.pl -lib xml.porphyromonas_gingivalis_w83.001.gz
 tracedb-to-frg.pl -frg xml.porphyromonas_gingivalis_w83.001.gz
 =>
 porphyromonas_gingivalis_w83.1.lib.frg  
 porphyromonas_gingivalis_w83.2.001.frg.bz2  
 porphyromonas_gingivalis_w83.3.lkg.frg 

• One lib at a time

 tracedb-to-frg.pl -only T13146 -xml xml.porphyromonas_gingivalis_w83.001.gz >& ! tracedb-to-frg.log
 tracedb-to-frg.pl -only T13146 -lib xml.porphyromonas_gingivalis_w83.001.gz >>& tracedb-to-frg.log
 tracedb-to-frg.pl -only T13146 -frg xml.porphyromonas_gingivalis_w83.001.gz >>& tracedb-to-frg.log

Illumina formatting

• use fastqToCA (wgs-6.0)

 fastqToCA \
   -insertsize 3000 300 \
   -libraryname 1 \
   -type illumina \
   -fastq /fs/szattic-asmg4/Bees/Bombus_impatiens/s_1_1_sequence.txt,/fs/szattic-asmg4/Bees/Bombus_impatiens/s_1_2_sequence.txt | \
    sed 's/ori:I/ori:O/' \
   > 1.frg
 =>
 {LIB                                                                                                                                      
 act:A                                                                                                                                     
 acc:1                                                                                                                                     
 ori:O                                                                                                                                     
 mea:3000.000                                                                                                                              
 std:300.000                                                                                                                               
 src:                                                                                                                                      
 .                                                                                                                                         
 nft:10                                                                                                                                    
 fea:                                                                                                                                      
 forceBOGunitigger=1                                                                                                                       
 isNotRandom=0                                                                                                                             
 doNotTrustHomopolymerRuns=0                                                                                                               
 doRemoveDuplicateReads=0                                                                                                                  
 doNotQVTrim=0                                                                                                                             
 goodBadQVThreshold=0                                                                                                                      
 doNotOverlapTrim=0                                                                                                                        
 usePackedFragments=1                                                                                                                      
 illuminaFastQType=illumina                                                                                                                
 illuminaSequence=/fs/szattic-asmg4/Bees/Bombus_impatiens/s_1_1_sequence.txt,/fs/szattic-asmg4/Bees/Bombus_impatiens/s_1_2_sequence.txt    
 .                                                                                                                                         
 }

 gatekeeper \
   -o asm.gkpStore.BUILDING \
   -T -F  \
   -E asm.gkpStore.errorLog \
   1.frg 
 => Segmentation fault

FASTQ

• This won't remove the linker:

 zcat SRR001351.fastq.gz | ~/bin/fastq2seq.pl  > SRR001351.seq
 zcat SRR001351.fastq.gz | ~/bin/fastq2qual.pl > SRR001351.qual
 convert-fasta-to-v2.pl -454 -l SRR001351 -s SRR001351.seq -q SRR001351.qual

• Cleaning seq

 ~/bin/OLD/cleanFastaq.pl \
    s_2_1_1.1kb_sequence.txt s_2_2_1.1kb_sequence.txt \                                     # input
    s_2_1_1.1kb_sequence.clean.txt  s_2_2_1.1kb_sequence.clean.txt \                        # output mates
    s_2_1_1.1kb_sequence.clean.single.txt  s_2_2_1.1kb_sequence.clean.single.txt            # output single

   echo prefix_1.fastq prefix_2.fastq  > prefix.ls
   quake.py -f prefix.ls

• Illumina note:

 We recommend removing reads that do not pass the GA analysis software Failed_Chastity fi lter before attempting to assemble the sequence.
 The chastity of a base call is the ratio of the intensity of the greatest signal divided by the sum of the two greatest signals. 
 Reads do not pass the quality fi lter if there are two or more base calls with chastity of less than 0.6 in the fi rst 25 cycles. 
 These reads have an “N” in the last column of the GA analysis software export file.

 Removed all reads that contained ambiguous characters (Ns).
 Removed reads that did not contain at least 25 Q30 bases among the first 35 cycles (s35). 

SFF

 sffToCA -libraryname SRR001351 -output SRR001351.frg SRR001351.sff

Running

• Command

 ~/wgs/Linux-amd64/bin/runCA -p pging -d testassembly porphyromonas_gingivalis_w83.* >& runCA.log

CAUG notes

Illimina mates

• short: innies PE (paired ends)
• long: outies MP (mate pairs) "lots of issues" (eg> : linker is one bp; solution : take 36bp from each end of the read to amke sure you don't get the linker

• wgs 6.1 : Illumina limit of 104 bp; can recompile the code o make it 232

fastqToCA -insertsize 280 30 -libraryname yersinia_pestis.100bp.0280bp.005x -type illumina -innie -fastq $PWD/yersinia_pestis.100bp.0280bp.005x.1.fastq,$PWD/yersinia_pestis.100bp.0280bp.005x.2.fastq > illumina.frg
{VER
ver:2
}
{LIB
act:A
acc:yersinia_pestis.100bp.0280bp.005x
ori:I
mea:280.000
std:30.000
src:
.
nft:12
fea:
forceBOGunitigger=1
shortOverlapModel=1
isNotRandom=0
doNotTrustHomopolymerRuns=0
doRemoveDuplicateReads=0
doNotQVTrim=0
goodBadQVThreshold=0
doNotOverlapTrim=0
usePackedFragments=1
illuminaFastQType=illumina
illuminaOrientation=innie
illuminaSequence=/home/training/CAUG/USERS/Daniela_Puiu/yersinia_pestis.100bp.0280bp.005x.1.fastq,/home/training/CAUG/USERS/Daniela_Puiu/yersinia_pestis.100bp.0280bp.005x.2.fastq
.
}
{VER
ver:1
}

#####################

gatekeeper -o illumina.gkpStore illumina.frg 
Starting file 'illumina.frg' at line 0.
Type set to ILLUMINA 1.3+.
Orientation set to INNIE.
Processing illumina reads from '/home/training/CAUG/USERS/Daniela_Puiu/yersinia_pestis.100bp.0280bp.005x.1.fastq'
                           and '/home/training/CAUG/USERS/Daniela_Puiu/yersinia_pestis.100bp.0280bp.005x.2.fastq'
GKP finished with no errors.

#####################

gatekeeper -dumpinfo illumina.gkpStore/
LOAD STATS

1       libInput
1       libLoaded
0       libErrors
0       libWarnings

0       frgInput
228742  frgLoaded
0       frgErrors
0       frgWarnings

0       lkgInput
0       lkgLoaded
0       lkgErrors
0       lkgWarnings

0       sffInput
0       sffLoaded
0       sffErrors
0       sffWarnings

0       sffLibCreated

0       plcInput
0       plcLoaded
0       plcErrors
0       plcWarnings

228742  numRandom
228742  numPacked
0       numNormal
0       numStrobe

GLOBAL STATS

228742  FRG
1       LIB

LibraryName                             numActiveFRG    numDeletedFRG   numMatedFRG     readLength      clearLength
GLOBAL                                  228742          0               218724          22874200        22191085
LegacyUnmatedReads                      0               0               0               0               0
yersinia_pestis.100bp.0280bp.005x       228742          0               218724          22874200        22191085

################

sffToCA -insertsize 3200 900 -libraryname porphyromonas_gingivalis_w83.3200bp.0900bp.E8YURXS01.1x.sff -trim chop -linker flx -output 454 porphyromonas_gingivalis_w83.3200bp.0900bp.E8YURXS01.1x.sff 
loadSFF()-- Loading 'porphyromonas_gingivalis_w83.3200bp.0900bp.E8YURXS01.1x.sff'.
removeDuplicateReads()-- from 1 to 19622
removeDuplicateReads()-- finished
detectMates()-- from 1 to 19621

#################

Overlap: check the hash load 
  6% is too low
~75% would be ideal