Bos taurus 3.0: Difference between revisions
Jump to navigation
Jump to search
Line 365: | Line 365: | ||
* nucmer params: -l 12 -c 65 -g 1000 -b 1000 | * nucmer params: -l 12 -c 65 -g 1000 -b 1000 | ||
* 24 * 30 = 720 alignments | * 24 * 30 = 720 alignments (except for BtChrU) | ||
* Best alignment summary: | * '''Alignments counts''' | ||
HsChr-BtChr.delta 532,866 | |||
HsChr-BtChr.filter-1.delta 392,789 | |||
HsChr-BtChr.filter-2K.delta 39,663 | |||
HsChr-BtChr.filter-2K-1.delta 38,185 | |||
HsChr-BtChr.filter-5K.delta 3,570 | |||
HsChr-BtChr.filter-5K-1.delta 3,560 | |||
* '''Best alignment summary:''' | |||
43 chr sets align over 1% | |||
nl #HsChr BtChr HsChrLen BtChrLen #al #2Kal min max median sum HsChrAl% BtChAl% | nl #HsChr BtChr HsChrLen BtChrLen #al #2Kal min max median sum HsChrAl% BtChAl% | ||
1 1 2 247249719 137060424 3229 238 13 9112 618 2763498 1 2 | 1 1 2 247249719 137060424 3229 238 13 9112 618 2763498 1 2 | ||
Line 412: | Line 422: | ||
42 22 17 49691432 75158596 1681 114 15 14130 558 1328105 2 1 | 42 22 17 49691432 75158596 1681 114 15 14130 558 1328105 2 1 | ||
43 23 30 154913754 148823899 17707 2895 15 19690 934 22575384 14 15 | 43 23 30 154913754 148823899 17707 2895 15 19690 934 22575384 14 15 | ||
... | |||
. . . . . 532794 41843 12 34598 718 493618379 . . | . . . . . 532794 41843 12 34598 718 493618379 . . | ||
* 54 chr sets have at least one 5K alignments | * 54 chr sets have at least one 5K alignments |
Revision as of 14:39, 14 September 2009
Sequence
The genome of the domestic cow, Bos taurus, was sequenced using a mixture of hierarchical and whole-genome shotgun sequencing methods.
Read download
- All reads were downloaded from the NCBI Trace Archive (TA) ftp: ftp://ftp.ncbi.nih.gov/pub/TraceDB/bos_taurus/
- There were 37,829,394 reads organized into 91 volumes
- 36,820,485 WGS, SHOTGUN, CLONEEND & FINISHING reads
- 36,170,352 quality reads
- 650,133 quality-less reads
- 1,008,909 EST & PCR reads
- 36,820,485 WGS, SHOTGUN, CLONEEND & FINISHING reads
- 25,312 read libraries
Sequencing centers
- Most reads were sequenced by the Baylor College of Medicine
TRACE_COUNT CENTER_NAME 1 35629020 BCM Baylor College of Medicine 2 737900 NISC NIH Intramural Sequencing Center 3 652614 BCCAGSC British Columbia Cancer Agency Genome Sciences Center 4 378871 MARC USDA, ARS, US Meat Animal Research Center 5 114753 UIUC University of Illinois at Urbana-Champaign 6 107367 BARC USDA, ARS, Beltsville Agricultural Research Center 7 65171 TIGR The Institute for Genome Research 8 53556 GSC Genoscope 9 43033 CENARGEN Embrapa Genetic Resources and Biotechnology 10 18623 SC The Sanger Center 11 15301 UOKNOR University of Oklahoma Norman Campus, Advanced Center for Genome Technology 12 10651 TIGR_JCVIJTC The Institute for Genomic Research, Traces generated at JCVIJTC 13 2485 UIACBCB University of Iowa Center for Bioinformatics and Computation Biology (UIACBCB) 14 49 WUGSC Washington University, Genome Sequencing Center 37829394 total total
Trace counts
TRACE_COUNT CENTER_NAME TRACE_TYPE_CODE 1 24863599 BCM* WGS 2 10748529 BCM* SHOTGUN 3 737900 NISC SHOTGUN 4 125597 BCCAGSC CLONEEND 5 114753 UIUC CLONEEND 6 65171 TIGR CLONEEND 7 53556 GSC CLONEEND 8 26246 CENARGEN WGS 9 25454 BARC CLONEEND 10 16892 BCM* CLONEEND 11 16787 CENARGEN CLONEEND 12 15150 UOKNOR SHOTGUN 13 10651 TIGR_JCVIJTC CLONEEND 14 151 UOKNOR FINISHING 15 49 WUGSC CLONEEND 36820485 total 16 527017 BCCAGSC EST 17 207204 MARC EST 18 171667 MARC PCR 19 81913 BARC EST 20 18623 SC EST 21 2485 UIACBCB EST 1008909 total
Data processing
Data issues
Issues:
- Qualities
- 650,133 reads don't have quality values and can't be reliably trimmed
- Libraries
- There are totally 25,312 libraries
- Very fragmented especially the SHOTGUN and CLONEEND ones; can't be accurately re-estimated by the assembler
- Clear ranges
- Many traces are missing vector trimming coordinates (CLV=CLIP_VECTOR_LEFT..CLIP_VECTOR_RIGHT) or don't contain 3' trimming information (CLIP_VECTOR_RIGHT==0)
- The read CLV's are need by the Celera Assembler overlap based trimming module (OBT) as input
- Solution: identify the sequencing vector & linker sequences for each library and re-trim the reads
Identify linkers
For each library identify linker sequences:
- Separate forward/reverse reads
- Identify most frequent kmers (8mers,24mers)
- Check if kmers a overrepresented
- Verify if the most frequent 8mer is present in the top 10 most frequent 24mers
- Align 24mers (extend them by a few bp) => linker
Identify vectors
For each library identify vector sequences:
- Align linkers to the opposite strand sequences (nucmer -l 12 -c 24 -r)
- Extract the subsequences following to linker (50..150bp)
- Align the subsequences; if they align we've probably identified the vector
- Identify the vector name/id by alignment to the UniVec database (nucmer -l 12 -c 24)
- Check if the forward/reverse vector(s) are the same : we should find a common vector sequence; the UniVec alignments should be adjacent
- create the Lucy vector & splice files that contain the linker+vector sequences
Trimming
- Run Lucy on quality reads
- Get CLV statistics: depending on the library, the Lucy CLV is 20bp+ shorter than the original CLV
- Trim reads according to Lucy output CLV
- Align Lucy trimmed reads to linker,vector,splice site & UniVec (there should be no alignments)
- Method worked on BCM & NISC libraries (~ 98% of the reads)
- For the other reads use the factory clipping points
BCM reads
- linker:
>J01636.linker.fwd 27bp TCGAGTTCGACTGCAAGTAGTTCATCA >J01636.linker.rev 27bp CTAATCAGATGGTACAGTAGTTCATCA
- vector: J01636 E.coli lactose operon with lacI, lacZ, lacY and lacA genes (7477 bp)
- avg(original CLV) - avg(Lucy CLV)> 20bp (1015 vs 973 in quality WGS reads , ...)
NISC reads
- linker:
>NGB00080.linker.fwd 24bp TATCATCGCCACTGTGGTGGAATT >NGB00080.linker.rev 26bp GCTGAAGCTCCATGTGGTGGAATTCC
- vector NGB00080 (pOTW13 with linkers)
- avg(original CLV) - avg(Lucy CLV)> 20bp (771 vs 747)
Preliminary assembly
- Assembly version: wgs-5.2
- Use only quality reads
- Set read CLV to Lucy CLV or original CLV
- Set non random flag = 1 on all reads except for WGS ones
- Set obtMerThreshold = 200 (default 1000)
- Set doOBT = 1
Input
Reads=35,348,776 # WGS, SHOTGUN, CLONEEND & FINISHING quality reads Libraries=25,312 # mostly SHOTGUN and BARC.CLONEEND
Output
TotalScaffolds=66,141 MaxBasesInScaffolds=26,048,998 MeanBasesInScaffolds=40,861 TotalContigsInScaffolds=120,461 MaxContigLength=627,911 MeanContigLength=22,436 TotalDegenContigs=269,031 MaxDegenContig=33,824 SingletonReads=3,721,123
DeletedReads=421,379 (too short or zero CLR)
Preliminary assembly processing
Read clear ranges
- Quality reads: extract OBT CLR from gatekeeper store
- Qualityless reads:
- Align them to contigs (no degenerates) : nucmer -l 50 -c 200 -b 10 -g 5 -d 0.05
- Set CLR to the maximum alignment coordinates or 50..min(len,600)
- Reduce CLR if there are multiple N's or low complexity regions in the read
Contaminant search
Databases:
- Ecoli : 22 completed genomes + plasmids
- UniVec_Core 1,348 sequences : mostly cloning vectors & primers, avg 250bp long
- OtherVec: 100 other vector sequences (mostly complete), identified by aligning UMD2.0 contaminants to GenBank
- bos_taurus UMD2.0 contaminant : 4,813 whole contigs and 30,329 partial contigs identified by NCBI as contamination in UMD2.0; many partial contigs contained cow sequences as well
- Databases FASTA files:
/nfshomes/dpuiu/db/Ecoli.all /nfshomes/dpuiu/db/UniVec_Core /nfshomes/dpuiu/db/OtherVec /nfshomes/dpuiu/db/bos_taurus.UMD2.contaminant.fasta
Alignment parameters:
nucmer -maxmatch -l 40 -c 100 -b 10 -g 5 -d 0.05
Contig/degenerate counts:
- 2,951/1,266 aligned to Ecoli
- 5,387/1,908 aligned to UniVec_Core
- 5,657/1,963 aligned to OtherVec
Read/mate counts: TO BE DELETED
- 40,699/22,607 in contaminated regions
Library estimates
- Some library estimates are complete wrong
Example: BCM.SHOTGUN libraries listed as long (180Kbp mean) are all short (2-6Kbp mean)
- Extract library insert estimates; merge libraries sequenced by same center that have similar mean/std : 25,312 libs => 344 libs
- Assign new library ids; assign average means & stdevs to the libraries
Final assembly
- Assembly version: wgs-5.2
- Use all traces
- Set read CLR to:
- Quality reads: OBT CLR
- Qualityless reads: alignment coordinates or 50..min(len,600)
- Set nonRandom flag = 1 on all reads except for WGS reads
- Set deleted flag = 1 on all reads deleted by OBT in the preliminary assembly
- Set obtMerThreshold = 200 (default 1000)
- Set doOBT = 0 (reads have been already trimmed)
Input
Reads=35,973,728 # WGS, SHOTGUN, CLONEEND & FINISHING with and without qualities Libraries=344
Output
TotalScaffolds=39,978 TotalContigsInScaffolds=90,135 MeanBasesInScaffolds=66,947 MaxBasesInScaffolds=3,3907,885 TotalContigsInScaffolds=90,135 MeanContigLength=29,693 MaxContigLength=1,160,130 TotalDegenContigs=251,413 MaxDegenContig=39,964 SingletonReads=3,634,305(10.24%)
Final assembly processing
Contaminant search
- Use same databases and alignment parameters as in preliminary assembly processing
- Delete full contaminants & trim partial contaminants
Delete summary:
- 65 Acinetobacter ctgs
- 91 other contaminant ctgs <2Kbp
- Total: 156 ctgs, 152 scf, 4105 reads
Trim summary:
- 12 contigs >=2Kbp , 44 reads
Marker mapping
- 126,013 total markers
- Avg distance between markers is 25Kbp; marker position error is 50Kbp
- Markers were aligned to all contigs/degenerates
- Best alignments with %IDY>90 & %Matched>85 were identified
- 107,271 markers align to 31,407 ctg & 2,640 scf
- 552 scf have markers from multiple chromosomes
- 212 scf have multiple markers from multiple chromosomes
- 38 scf have multiple adjacent markers from multiple chromosomes: MIGHT BE MISASSEMBLED
- 628 markers align to 562 degenerates
Scaffold/contig breaking
- Analyze 38 scf that have multiple adjacent markers from multiple chromosomes
- Compute coverage in the suspicious region (between different chromosome markers):
- read cvg
- mate ctg: good, bad
- Break ctg/scf unless the region has "high read cvg" , "high good mate cvg" , "low bad mate cvg"
- Break summary:
- 14 scaffolds
- 15 breaks : 8 on the same contig , 3 on adjacent contigs , 4 on non adjacent contigs
Assignment to chromosomes
Markers
- 2640 scaffolds and 562 degenerates have markers
- Assignment to chromosomes: use best alignment & majority rule
- Position:
- Filter out outliers according to position on chromosome & scaffold (interquartile range method)
- Compute the average position on chromosome of the markers
- Orientation:
- use LeastSequareFit method : if slope is positive => forward; if slope is negative => reverse
- if only 1 markers/scaffolds => direction=unknown (0)
Human synteny
- Align all scaffolds/degenerates to the 24 Human chromosomes; filter all alignments longer than 200bp
nucmer -mum -l 12 -c 30 -g 1000 delta-filter -q -l 200
- 9,914 scaffolds and 16,527 degenerates align to Human chromosomes; most alignments are short, just over 200bp
Combine Human synteny & Marker data
- 1,908 scaffolds and 118 degenerates both align to human and contain markers
- 10,790 scaffolds and 16,590 degenerates align to human or contain markers
- Try to infer the position/orientation on the chromosomes for the scaffolds/degenerates that align to human but contain no markers
- Iteratively:
- Find 2 adjacent scaffolds (preferably on left & right side) which both align to human, contain markers and placements agree (chromosome, position, direction)
- Otherwise, find 1 adjacent scaffolds which both aligns to human, and contains markers
- Extrapolate the position/orientation of the "unplaced" sequence based on its neighbor(s)
- Sort the scaffolds/degenerates based on chromosome positions, identify incorrect markers & alignments, remove them from the input data and repeat the process
By linking information
- Once scaffolds/degenerates were assign to chromosome use mate pair information to refine placements
- Identify unplaced scaffolds/degenerates linked to placed scaffolds/degenerates and fit them into gaps
Comparison to UMD2.0
Alignment parameters:
nucmer -mum -l 200 -c 1000
Haplotype search
Daniela:
- Place scf/deg on Chr
- Align each pair X,Y (len(X)<len(Y))of adjacent/overlapping scf/deg : nucmer -mum -l 40 -c 250 ( => avg 96 %id)
- compute (X,Y) cvg
- identify the X regions which had no alignments to Y; if the length of these regions were less than 2K bp => X is a variant
Guillaume:
- there is a contig Y larger than X and, X and Y are placed on top of each other (with some play allowed)
- there are a high quality sequence alignment such that: at least 200 bases out of the first 400 bases of X align with Y AND at least 200 bases out of the last 400 bases of X align with Y.
- In other word, the ends of X have to align well with Y, but the middle can be significantly different.
Files:
/fs/szasmg3/bos_taurus/UMD_Freeze3.0/contigs.haplotype-variants.fa.gz # 40611 haplotype-variants sequences /fs/szasmg3/bos_taurus/UMD_Freeze3.0/contigs.haplotype-variants.ids # 40611 haplotype-variants sequence ids /fs/szasmg3/bos_taurus/UMD_Freeze3.0/contigs.haplotype-variants.pairs # 40300 pairs (haplotype-variants & reference sequences) /fs/szasmg3/bos_taurus/UMD_Freeze3.0/contigs.haplotype-variants.pairs.delta # 39665 alignment pairs (haplotype-variant is the query) /fs/szasmg3/bos_taurus/UMD_Freeze3.0/contigs.haplotype-variants.pairs.cvg # 39665 coverage pairs (ref=col 1 ; haplotype-variant=col 5)
Summary:
. elem <=0 >0 min max mean med n50 sum ctg+deg 40611 0 40611 263 97877 1476 1205 1372 59958728 ctg 29452 0 29452 471 97877 1631 1297 1469 48039280 deg 11159 0 11159 263 12208 1068 979 1006 11919448
Other Files
/scratch1/bos_taurus/Assembly/2009_0312_CA/scf_placements/UMD_Freeze3.0/ 39864 contigs.haplotype-variants.daniela.pairs 443 contigs.haplotype-variants.guillaume.pairs 436 contigs.haplotype-variants.guillaume.pairs.orig 334 contigs.haplotype-variants.ids.missing
Issues:
. elem <=0 >0 min max mean med n50 sum missing 339 0 339 462 97877 1613 1011 1204 547135 mislabeled 6 0 6 2973 97877 31271 8275 97877 187628
Mislabeled haplotypes:
Chr begin end Pos W ctg 1 len dir Chr2 131428091 131475109 5387 W 7180001925346 1 47019 + Chr3 11270384 11278308 637 W 7180002020315 1 7925 - Chr12 9183395 9281271 429 W 7180002024890 1 97877 + Chr14 34404395 34412669 3579 W 7180002021388 1 8275 - Chr15 50865607 50889165 2531 W 7180002015261 1 23559 - ChrU 8589989 8592961 5791 W 7180002026074 1 2973 +
Chromosome mapping
Assembly Summary
...
Hs vs Bt
- nucmer params: -l 12 -c 65 -g 1000 -b 1000
- 24 * 30 = 720 alignments (except for BtChrU)
- Alignments counts
HsChr-BtChr.delta 532,866 HsChr-BtChr.filter-1.delta 392,789 HsChr-BtChr.filter-2K.delta 39,663 HsChr-BtChr.filter-2K-1.delta 38,185 HsChr-BtChr.filter-5K.delta 3,570 HsChr-BtChr.filter-5K-1.delta 3,560
- Best alignment summary:
43 chr sets align over 1% nl #HsChr BtChr HsChrLen BtChrLen #al #2Kal min max median sum HsChrAl% BtChAl% 1 1 2 247249719 137060424 3229 238 13 9112 618 2763498 1 2 2 1 3 247249719 121430405 18208 2573 15 19711 855 21390308 8 17 3 1 16 247249719 81724687 11817 1057 13 26957 727 11384863 4 13 4 2 2 242951149 137060424 18147 2379 15 28727 857 21057071 8 15 5 2 11 242951149 107310763 13922 1292 15 21608 740 13879256 5 12 6 3 1 199501827 158337067 17808 2100 14 24243 837 19561651 9 12 7 3 22 199501827 61435874 10408 973 14 18031 752 10518104 5 17 8 4 6 191273063 119458736 14857 1378 14 15997 776 14933999 7 12 9 4 17 191273063 75158596 5498 507 15 14456 800 5603523 2 7 10 5 7 180857866 112638659 14201 1891 13 26095 866 16433557 9 14 11 5 20 180857866 72042655 9426 855 15 13811 767 9274385 5 12 12 6 9 170899992 105708250 14577 1559 14 18079 818 15402674 9 14 13 6 23 170899992 52530062 8082 804 14 15593 731 8040414 4 15 14 7 4 158821424 120829699 17230 1835 13 34598 803 18223663 11 15 15 7 25 158821424 42904170 2574 129 14 11489 536 1927855 1 4 16 8 8 146274826 113384836 1891 100 15 11015 598 1491998 1 1 17 8 14 146274826 84648390 11817 1025 14 30396 744 11637718 7 13 18 8 27 146274826 45407902 3207 143 14 8051 651 2597335 1 5 19 9 8 140273252 113384836 14511 1540 14 17317 801 15339451 10 13 20 9 11 140273252 107310763 2874 253 15 19675 609 2568633 1 2 21 10 13 135374737 84240350 4371 245 14 12436 691 3784917 2 4 22 10 26 135374737 51681464 8467 917 13 19568 757 8745051 6 16 23 10 28 135374737 46312546 5712 567 14 11567 761 5753856 4 12 24 11 15 134452384 85296676 12655 1340 14 22774 808 13419432 9 15 25 11 29 134452384 51505224 7120 528 13 16632 685 6446478 4 12 26 12 5 132349534 121191424 15222 1787 14 29685 804 16677907 12 13 27 12 17 132349534 75158596 3924 213 13 8484 585 3042898 2 4 28 13 12 114142980 91163125 11389 837 14 19959 738 10665467 9 11 29 14 10 106368585 104305016 8626 980 14 17863 807 9394454 8 9 30 14 21 106368585 71599096 5319 651 14 18953 785 5766996 5 8 31 15 10 100338915 104305016 7914 1088 14 16238 805 8977632 8 8 32 15 21 100338915 71599096 5057 438 15 10071 686 4689441 4 6 33 16 18 88827254 66004023 7632 655 14 19289 711 7291016 8 11 34 16 25 88827254 42904170 4680 278 12 9947 589 3698904 4 8 35 17 19 78774742 64057457 11876 1195 13 15954 651 11528275 14 17 36 18 24 76117153 62714930 9815 755 13 16525 704 9096937 11 14 37 19 7 63811651 112638659 2718 181 14 9082 507 2027086 3 1 38 19 18 63811651 66004023 3765 354 12 13895 624 3457684 5 5 39 20 13 62435964 84240350 8838 692 13 13905 671 7962468 12 9 40 21 1 46944323 158337067 3653 205 13 10470 644 3032218 6 1 41 22 5 49691432 121191424 2618 161 14 7411 540 1966904 3 1 42 22 17 49691432 75158596 1681 114 15 14130 558 1328105 2 1 43 23 30 154913754 148823899 17707 2895 15 19690 934 22575384 14 15 ... . . . . . 532794 41843 12 34598 718 493618379 . .
- 54 chr sets have at least one 5K alignments