This is a supplementary page of Genome-wide mutagenesis of Zea mays L. using RescueMu. Here you will be able to find the detailed materials and methods that used to assemble and analyze RescueMu GSS.

• Materials
• Assembly
  • Procedure
  • Result
• Analysis
  • Genic regions
  • Repeat regions
  • Genetic mapping
  • Insertion site profile

Materials

• RescueMuGSSraw.fsa.gz: Raw RescueMu GSS provided from Dr. Walbot's lab at Stanford.
• RescueMuGSS.fsa.gz: RescueMu GSS downloaded from GenBank (search GenBank sequences with "Zea mays [ORGN] AND GSS [KYWD] AND Walbot [AUTH]").

When RescueMu inserts into maize genome, it usually creates a 9 bp host target site duplication (TSD) next to its terminal inverted repeats (TIRs) at each end of the transposon. After plasmid rescue, sequencing was performed using primers complementary to the left or right TIRs. To re-construct the original genome sequence at a particular RescueMu insertion site, we need to identify the 9 bp TSDs and use them to connect the left and right sequences flanking a RescueMu insertion site. For sequences derived from the same sequencing plasmid, the reverse-complementary of the left should overlap the right sequence starting with the 9 bp TSD. The raw RescueMu GSS contains some TIR sequences, but these were masked with Xs using CrossMatch program. In addition, the raw sequences contain low quality regions (phred quality score < 20) that were trimmed off prior to submission to GenBank.

• What is the total number of Raw sequences: 98,111 (How can you get this number? Use linux command "grep -c '>' RescueMuGSSraw.fsa")
• What is the total number of RescueMu plasmid in Raw sequences? 57,140 (How can you get this number? Use linux command " grep '>' RescueMuGSSraw.fsa | awk '{print $1}' | cut -d '.' -f1 | sort -u | wc -l")
• What is the total number of GenBank sequences? 131,364 (How can you get this number? Use linux command "grep -c '>' RescueMuGSS.fsa")
• What is the total number of RescueMu plasmid in GenBank? 57,022 (How can you get this number? Use linux command "grep '>' RescueMuGSS.fsa | cut -d ' ' -f2 | cut -d '.' -f 1 | sort -u | wc -l"). Note that there are data from numerous additional grids that were partially sequenced or for which sequencing is still in progress; these grids were not included in the analysis, but the same tools we developed to analyze RescueMu for the publication can be used with the additional sequences.

Why is the number of GenBank sequences larger than the raw sequences count? This occurs because one raw sequence (corresponding to one sequencing read) is sometimes split into two sequences for separate submission into GenBank. This procedure was used to identify sequences that might be made contiguous during the plasmid rescue circularization of a segment of maize genomic DNA. For example, the sequence 1006001B06.x1 in the raw sequence was split into two sequences at GenBank (1006001B06.x1 and 1006001B06.1EL_x1). The split occurs if either BamHI, BglII, or their possible mix ligation site was encountered. More information about the splitting can be found at http://www.mutransposon.org/project/RescueMu/zmdb/library-plate/GridGprep.html. As you can also see, a minor number of RescueMu plasmids didn't submit into GenBank due to low qualities.

• In GenBank, what is the total number of X sequences (corresponding to the right sequences)? 49,577 (How can you get this number? Use linux command "grep -c '\.x' RescueMuGSS.fsa")
• In GenBank, what is the total number of ELX sequences (corresponding to the split of right sequences)? 23,023 (How can you get this number? Use linux command "grep -c 'EL_x' RescueMuGSS.fsa")
• In GenBank, what is the total number of Y sequences (corresponding to the left sequences): 40,887 (How can you get this number? Use linux command "grep -c '\.y' RescueMuGSS.fsa"). Note that the primer used for Y sequencing was less efficient than the right primer, therefore there are fewer sequences.
• In GenBank, what is the total number of ELY sequences (corresponding to the split of left sequences): 17,877 (How can you get this number? Use linux command "grep -c 'EL_y' RescueMuGSS.fsa")

Assembly

Step 2 - repeat sequence masking: The vast majority of the maize genome was estimated to be repeat elements (e.g., retrotransposons). Our stragety here is to mask the repeat elements to avoid different loci from being clustered together due to their common repeat elements. Therefore, we screened the above 130,861 vector-free sequences against maize repeat elements annotated by TIGR. The screening was performed by the VMATCH program (-d -p -l 50 -h 3 -identity 95). The screened repeat sequences were trimmed off at the ends by a perl script (RepeatTrim.pl) and partitioned into the following two pools:

• The major one: 127,708 repeat-free sequences - this data set represents the potential "userful" sequences from the RescueMu approach. Our assembly is focused on this data set in order to evaluate the gene discovery efficiency.
• The minor one: 3,153 repeat-containing sequences - this data set represents the putative "junk" sequences from the RescueMu approach. We still assembled this data set in order to estimate the number of repeat regions RescueMu inserted into. however, the number of resulted contigs is unreliable and likely an underestimate of the true number recovered because of the intrinsic difficulty with assembling repetitive DNA.

Step 3 - identify parental sequences: Any given RescueMu-transformed plant contain the parental RescueMu elements that were recovered at a high frequency during sequencing (from every sequenced row or column). Because our goal is to analyze the gene discovery by newly inserted RescueMu (i.e., where the non-parentals inserted into the maize genome), we decide to filter out the parental sequences as much as possible. We designed a fast-filter to quickly identify large duplicated (or near-identical) sequences for each grid based on the property that parental sequences were recovered from every rows/columns of grid plants. The following is how our filter works:

1. From the above 127,708 repeat-free sequences, extract the left (Y) and right (X) sequences from each grid separately by a perl script (ExtractGridXYseq.pl).
   
   Sequence count from the above 127,708 repeat-free sequences:
   • Total number of X sequences from Grid G: 10,340 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid G' | wc -l)
   • Total number of Y sequences from Grid G: 8,000 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid G' | wc -l)
   • Total number of X sequences from Grid H: 10,059 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid H' | wc -l)
   • Total number of Y sequences from Grid H: 8,450 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid H' | wc -l)
   • Total number of X sequences from Grid I: 9,427 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid I' | wc -l)
   • Total number of Y sequences from Grid I: 7,764 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid I' | wc -l)
   • Total number of X sequences from Grid K: 4,851 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid K' | wc -l)
   • Total number of Y sequences from Grid K: 4,348 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid K' | wc -l)
   • Total number of X sequences from Grid M: 3,402 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid M' | wc -l)
   • Total number of Y sequences from Grid M: 3,074 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid M' | wc -l)
   • Total number of X sequences from Grid P: 10,628 (grep '\.x' TIGRrepeatTrim_RescueMuGSS | grep 'Grid P' | wc -l)
   • Total number of Y sequences from Grid p: 9,075 (grep '\.y' TIGRrepeatTrim_RescueMuGSS | grep 'Grid P' | wc -l)
   
   
   
2. Single-linkage cluster the above X, Y sequences from each grid (GridG_XYseq.fsa.gz, GridH_XYseq.fsa.gz, GridI_XYseq.fsa.gz, GridK_XYseq.fsa.gz, GridM_XYseq.fsa.gz, GridP_XYseq.fsa.gz) by VMATCH program (-d -p -v -l 100 -seedlength 50 -exdrop 1 -identity 99 -dbcluster 98 95). The clustering was performed as the following: given two sequences, first checking if they share 50 bp identical regions (seeds); if so, extend that seeds to both directions, allowing only one mismatch/gap during the extension. Then the two sequences would be considered as "neighbors" if they satisfy the following three criteria: (1) their match region is no less than 100 bp, (2) >=99% identity in the match region; (3) the match region is 98% of the total length of the shorter sequence, 95% of the longer sequence. The criteria clusters together extremely similar (or identical) sequences together. Single-linkage cluster: a member was added into a cluster if this member is a neighbor to at least one other member in the cluster.



3. Then a perl script (ParentalVmatchCluster.pl) was used to identify parental VMATCH clusters.
   
   

   Table 1: Grid Sequencing Information

   Grid	Num of Sequenced Rows	Num of Sequenced Cols
   G	45	1
   H	1	35
   I	36	2
   K	27	3
   M	0	31
   P	42	1

   
   Some key points of the ParentalVmatchCluster.pl script: the above table (showing the number of sequenced rows/columns for each grid) and the RescueMu sequence grid location information were combined to identify parentals. If a VMATCH cluster contains plasmids that were derived from all the rows (or columns, whichever the one is larger), then this VMATCH cluster is marked as parental cluster in which all the plasmids were marked as parental plasmids. For example, if a grid G cluster has plasmids from 45 different rows, then this is a parental cluster at grid G. Grid P is a special case: the parentals from row 1 to row 37 are different from the rest rows; therefore we only counted clusters with plasmids from row 1 to row 37 as parentals.
   
   The result: total 16,122 plasmids were identified as parentals plasmids from all six grids, representing total 42,738 repeat-free sequences that were excluded from the subsequent assembly. Note that these 16,122 parental plasmids were identified by very stringent matching criteria. Therefore, some parental sequences may slip through this parental filter. Please read on to see how we calibrate that.

Step 4 - plasmid-level assembly: As stated above, each sequencing plasmid may generate up to four sequences that were separately submitted to the GenBank (X, Y, EL_X, EL_Y). Therefore, before we start assembling sequences from all the plasmids, we first did a plasmid-level assembly to (1) connect X and Y together through 9 bp TSD; (2) assign orientation of EL sequences to X or Y if they share enough sequence overlap (e.g., if EL_Y overlaps with X, it indicates EL_Y and X were on the same side of the RescueMu ligation). The plasmid assembly was carried out by a perl script, RescueMuNonParentalPlasmidAssembly.pl.

Some key points for this perl script:

• Duplicated sequencing: In some cases, a plasmid was sequenced more than once from the same primer, for example from the right sequencing primer to obtain longer sequence, all the resulted sequences (including the short ones) were labeled as .x1, .x2 etc and submitted to the GenBank. In such cases, only the longest sequence was chosen to participate this plasmid-level assembly. For example, if x1 is longer than x2, we only used x1 for the assembly.
• Determine the 9 bp TSD. For plasmids generating both the left and right sequences, 9 bp TSD can be only identified if (1) the TIRs in the raw sequences were clearly masked; (2) at least 8 out of 9 bases are identical for the left and right sequences immediately outside the TIRs. If the 9 bp from left and right do not match perfectly (i.e., only 8/9 are identical), we assign the mismatch position by looking into the corresponding GenBank sequences for sequence quality. For example, if the left has high quality base on that position while the right does not, we use the base from the left. If both have high quality bases, we arbitrarily choose one from either the left or right.
• Assemble EL sequence. If the 9 bp TSD can be determined to connect X and Y (XYtuc), then we try to match EL sequences to the XYtuc; otherwise, EL was matched to X or Y. The matching of EL was performed with CAP3 (a modified version of CAP3 by ourselves; -o 21 -p 95 -s 401).
• Sequence quality. We always use GenBank sequences (high quality; average phred 20 or fewer than one error in 3000 bases). The raw sequences were only used for identifying 9 bp TSD and connecting the left and right by the identified 9 bp host sequence direct repeat (TSD). The low quality region in the raw sequences are always being trimmed off to match the ones deposited in the GenBank, except the ones that around the 9 bp region because we need them to identify TSD). This means that in some cases, we had to use some low-quality sequences in the 9 bp regions to connect the left and right genomic regions flanking a specific RescueMu insertion site).
• Inconsistency. The GenBank sequences SHOULD be the substring of the corresponding Raw sequences. If not, we throw away such inconsistent sequences that are likely caused by tracking errors.
• Multiple EL. Sequences with an EL site (restriction site from BamHI, BglII or their mix ligation site; see http://www.zmdb.iastate.edu/zmdb/library-plate/GridGprep.html for detail explanation) presents the challenge of determining orientation with respect to the RescueMu insertion site. Sequences with multiple EL sites have greater uncertainty. To ensure our assembly quality, we decided to set aside such multiple EL sequences.

1. Sequences thrown away due to the inconsistency: 784 (0.6% of the total sequences to be assembled)
2. Sequences set aside due to multiple EL sites: 9,797 (7.5% of the total sequences to be assembled)
3. Xs and Ys from 8,327 plasmids were connected through 9 bp TSD.
4. The final result of this plasmid-level assembly is 66,905 putative Non-Parental Plasmid Assembly Sequences, derived from 38,526 plasmids.
   
   

• A series of perl scripts were developed to proceed through PaCE/CAP3 and process their results.
• A sample display of GSS contig - ZM_RM_GSStuc03-10-31.4765

Parentals re-visit: as stated above, we have identified total 16,122 parental plasmids and excluded them from the PaCE/CAP3 assembly. Here we want to re-visit this issue and answer the following two questions:

• How many insertion sites do those parental plasmids represent?
• How many parentals were not filtered and went into PaCE/CAP3 assembly?

1. Build representative contigs: For each grid, a perl script (RandomXYSeqCluster.pl was used to randomly select X and Y sequences from 50 different parental plasmids from each parental cluster. The X and Y were joined together through their 9 bp TSD. Then, all the selected sequences were pooled together for each grid and ran CAP3 (-p 90 -s 700) to generate contigs. For example, there were 2 parental clusters from grid G. Then 50 random sequences from each cluster were selected and combined (100 sequences) for CAP3 to generate a representative 1 contig. Notice that the parental cluster was built using individual X and Y sequences, which means that X and Y from the same plasmids were separated in different cluster because they only shared at most 9 bp TSD region. Now that we joined X and Y, those different cluster can be merged together (that is why we had 2 clusters from G merged into only 1 representative contig).
2. Extend representative contigs: We took the advantage of the existing huge amount of maize GSS sequences from the Methylation filtration and High Cot library to extend the contigs we built from (1).
   • Methylation filtration GSS contigs (MF.TUG.fasta) and High Cot GSS contig (HC.TUG.fasta).
   • Representative contigs from each grid.
   • Then the VMATCH program (-d -p -l 40 -h 2 -identity 95) was used to select any potential MF.TUG.fasta or HC.TUG.fasta that matches the representative contigs. The selected potentials were combined with the representative contigs to generate CAP3 contigs (-p 90 -s 700; extended contigs).
   
   
3. Verify parental sequences: The reason to build the above extended parental contigs is to verify all the parental sequences identified by our parental-filter. All the sequences that we claimed to be parental should be covered by such extended contigs. For each grid, we used BLASTN program (-e 1e-4) to compare the extended contigs with all the parental sequences from that grid.
   
4. Deep Clean Parentals: Now that we have the parental contigs. We can use those contigs to screen the above PaCE/CAP3 contigs (from the above Step 5) by using the BLASTN program (-e 1e-4). The result further identified 508 parental assemblies (Alignment; representing 16,331 GSS sequences from the PaCE/CAP3 assembly sequences). All the parental sequence Stanford IDs (42,738 from Step 3 + 16,331 from parentals deep clean; total 24,875 plasmids).

   Table 2: Parental plasmids

   Grid	Parental Plasmid Count
   G	2,048
   H	5,231
   I	5,246
   K	3,091
   M	1,422
   P	7,837

Repeats re-visit: as stated above, the GSSs were masked by maize repeat elements. The above contigs were assembled from repeat-trimmed GSS. In order to estimate the repeat contents from the original sequences obtained from RescueMu, we put back the trimmed repeats into the contigs. Specifically, 318 contigs contain repeat-trimmed GSS ans were re-assembled with CAP3 using the original GSS. The re-assembling resulted in 367 repeat-containing contigs. In addition, the above 3,153 repeat-containing GSSs (the minor pool) were also assembled into contigs.

Final Assembly results

• 25,068 final contig sequences with total size of 9.92 MB.
• The contigs were assembled from 62,803 GSSs derived from 32,639 plasmids.

Genomic loci: The above PaCE and CAP3 assembly is purely based on sequence overlap (i.e., if sequences share some significant overlap, they will be clustered and assembled into a contig). Then ideally, each assembly would correspond to a unique genomic location on the maize genome. Unfortunately, there are many cases where sequences that were indeed from the same genomic loci DO NOT get clustered together and are separated into different assembly. For example, contigs ZM_RM_GSStuc03-10-31.895 and ZM_RM_GSStuc03-10-31.896 do not share significant sequence overlap. But they both clearly represent the same genomic locus, because they were derived from sequences from the same plasmids. In this particular case, those sequences can not be clustered together because the 9 bp TSD can not be identified to connect the left (Y seq) and right (X seq).

Therefore, in addition to the sequence overlap, our analysis includes another single-linkage clustering using clone-pair constraints. Any contigs that share sequences from the same plasmids will be assigned to the same genomic locus. Such single-linkage clustering was carried out using a perl script. As a result, the above 25,068 contigs were clustered into 14,265 genomic loci (note: 13,187 loci were clustered from repeat-poor contigs and 1,078 loci were clustered from repeat-rich contigs). Single-linkage cluster: a member was added into a cluster if this member is a neighbor to at least one other member in the cluster.

Analysis

(A)

(B)

Genetic Mapping: All the ESTs that were genetically mapped (total 1,775 unique ESTs) on the maize genome were downloaded from the MaizeGDB. Their sequences were aligned to the above GSS contigs as described above. Total 131 matched GSS contigs were able to be plotted on IBM neighbor map coordinates system. The result clearly shows that RescueMu jumps to all the maize chromosomes as expected.

RescueMu Insertion Site: The RescueMu tagging approach generates not only maize genomic sequences but also potential maize mutants for functional analysis. We have located 3,999 independent RescueMu insertion sites on 3,688 GSS contigs. Note that in order to locate the RescueMu insertion site, the X and Y sequences from the same plasmids must be available and their 9 bp TSD could also be identified as described above. Since the vast majority of the sequencing plasmids do not have both X and Y sequences, we can not locate all the insertion site. However, these 3,686 GSS contigs contain 3,999 confirmed insertion site and still provide a good sample size for our estimation. Out of all the 3,999 insertion sites, 968 insertion sites are WITHIN confirmed gene structure (exon-intron structure) provided by EST/cDNA matching. 849 insertion sites are within exons (119 insertion sites are within introns), showing a strong perference for RescueMu to insert into exons, which is consistent with Mu being an efficient mutagen.