Computational protocol: Genome-wide characterization of intergenic polyadenylation sites redefines gene spaces in Arabidopsis thaliana

Similar protocols

Protocol publication

[…] Distances from each IPAC to its 5’ sense gene, 3’ sense gene, 5’ antisense gene and 3’ antisense gene were considered as the main metric to classify IPACs. The same number of positions in intergenic regions were randomly selected as background control to estimate FDR. Histograms (Additional file : Figures S2 and S3) were plotted to show the distributions of distances from IPACs or random intergenic positions to sense and antisense genes using an R package called ggplot2 [] with 100 breakpoints. To determine an appropriate distance for SE-IPACs, the numbers of IPACs and random background positions were calculated from the nearest annotated stop codon to a maximum distance of 2000 nt by an interval of 20 nt (Fig. ). FDR was defined as the ratio of background positions in a specific distance to the total number of background positions within the 2000 nt distance. SE-IPACs were first defined according to the distance from IPACs to their upstream genes. The distance that meets the criterion of FDR < 0.1 (700 nt) was obtained to differentiate SE-IPACs and other IPACs. IPACs that are located in their antisense genes were classified as A-IPACs. To further remove possible internally primed artifacts of A-IPACs, if the downstream 20 nt of a A-IPAC has more than 50 % of A or more than 60 % of A and G, then it is discarded. To determine the case of SO-IPACs, we first removed IPACs that are close to 3’ antisense genes using a similar strategy for SE-IPACs. When the distance to the stop codon is 200 nt, the FDR from background positions is 0.01. There are 791 IPACs located within the 200 nt region. These IPACs close to antisense promoter were removed. The rest IPACs of 5262 are either close to 3’ sense genes or far from all nearby genes. Setting the same distance cutoff as SE-IPACs, 344 IPACs are located within 700 nt from the corresponding 3’ sense gene. These IPACs close to sense promoter were also discarded for further analysis. The rest 1317 IPACs which are far from all nearby genes were classified as SO-IPACs. [...] EST data were used for the conserved estimation of novel transcripts identified from SO-IPACs. EST/cDNA data were downloaded from TAIR website (ftp://ftp.arabidopsis.org/home/tair/Sequences/). These sequences were mapped to TAIR10 reference genome using GMAP []. The coordinates of 3’-ends of mapped ESTs were recorded. 180,805 aligned ESTs were located in intergenic regions. To verify SO-IPACs, the coordinates of SO-IPACs and EST 3’-ends were compared. If a SO-IPAC is located within 500 nt of any aligned EST, then it is considered as EST-verified IPAC.Publicly available RNA-seq data [–] were also employed to validate SO-IPACs. For the data from Filichkin et al. [], reads sequenced from normal physiological conditions were used and the 76 nt reads were truncated to 32 nt as recommended in []. For the data from [, ], we directly used the TopHat alignment files available at www.compbio.dundee.ac.uk/polyADB/ []. We ran TopHat [] on each source of RNA-seq data using the following parameters: no multi-hits (−g 1), minimum anchor length 10 (−a 10), and minimum and maximum intron length 40 and 5000, respectively. Next, each alignment output from TopHat was assembled into putative transcripts by Cufflinks using TAIR10 genome annotation as the reference []. Then all Cufflinks assemblies were merged together using Cuffmerge [] and genes in the final assembly denoted by “CUFF” and not overlapped with any annotated gene model in TAIR 10 were considered as novel genes. If a SO-IPAC is located within the novel genes identified from the TopHat/Cufflinks pipeline, then it is validated by the RNA-seq data. Because both the SE-IPAC and A-IPAC were associated with known gene models, we used genes in the Cufflinks assembly denoted by “CUFF” with FKPM (fragments per kilobase of exon per million fragments mapped) >0 to verify these two classes of IPACs.Highly conserved blocks from 20 Angiosperm Plant Genomes [] were employed to study the conservation of IPACs across plant species. Conserved regions of the mostCons set were downloaded from the Arabidopsis genome browser (http://genome.genetics.rutgers.edu/). Coordinates of IPACs were compared with the start/end coordinates of the conserved regions to examine whether an IPAC is located in a conserved region or not.Gene functional annotation was performed using the DAVID Bioinformatics Resources []. For this, lists of genes associated with SE-IPACs and A-IPACs were compared with the list of all genes that are represented in the PAC datasets. Category designations: GOTERM_BP_FAT: Gene Ontology Biological Process classes; GOTERM_MF_FAT: Gene Ontology Molecular Function classes; GOTERM_CC_FAT: Gene Ontology Cell Component classes.To estimate the protein coding potential of SO-IPACs in our dataset, we applied the Coding Potential Calculator (CPC, []) scanning for evidence for protein coding capacity. The input is the dataset of sequences of upstream 500 nt to 200 nt of all SO-IPACs. The output is the CPC score for each input sequence. Wilcoxon rank sum test for 3’-UTR lengths between two groups was carried out by the function wilcox.test in program R (http://www.r-project.org). To test the relation between SO-IPAC and ORF or SSP, the predicted small ORFs were downloaded from http://peptidome.missouri.edu []. Because the coordinates of these ORFs were based on TAIR6 annotation, we remapped DNA sequences of the ORFs to the TAIR10 genome using Bowtie2 [] to obtain the updated coordinates. If a SO-IPAC is located within downstream 2000 nt of an ORF, then it is considered as associated with the ORF. To examine the correlation of expression level between A-IPAC and its target antisense gene, the total number of PATs of sense A-IPACs and targeted antisense genes were considered as expression levels. Test for correlation between paired samples based on Pearson's product moment correlation coefficient was performed by the function cor.test in program R. The occurrences of AATAAA and its variants in a studied region of IPACs were counted using a script written in Perl. […]

Pipeline specifications

Software tools Ggplot2, GMAP, TopHat, Cufflinks, DAVID, CPC, Bowtie2
Databases TAIR
Applications Miscellaneous, RNA-seq analysis, Transcription analysis, Genome data visualization
Organisms Arabidopsis thaliana