Bos taurus 3.0
Jump to navigation
Jump to search
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
. ctg+deg <2Kbp >=2Kbp min max mean med n50 sum ====================================================================================================== Chr1..29,X 72479 20862 51617 65 1160130 36424 12941 103785 2639984487 ChrU 3285 2404 881 224 179692 2890 1338 5425 9496583 Chr 75764 23266 52498 65 1160130 34970 11207 96955 2649481070 contigs.haplotype-variants 40611 36984 3627 263 97877 1476 1205 1372 59958728 deg.unplaced.less_2K 224933 224933 0 65 1996 972 983 990 218837572 ChrY-contigs 314 266 48 224 26490 2210 973 6539 694140 ChrY-contigs.SHOTGUN_ONLY 144 140 4 804 4224 993 882 888 143047 ======================================================================================================
Hs vs Bt
nucmer params:
-l 12 -c 65 -g 1000 -b 1000