class: center, middle, inverse, title-slide # Visualize thetas and Fst results ### Jinliang Yang ### Feb. 13th, 2020 --- # Visualize the results In local computer, using `R`: #### Barplot for SFS ```r s <- scan('../../cache/out.sfs') s <- s[-c(1,length(s))] s <- s/sum(s) barplot(s,names=1:length(s), main='SFS') ``` #### Histgram distribution of the theta values ```r hist(t$Pairwise) ``` #### Scatter plot of the Fst values ```r fst <- read.table("cache/fst_win.txt", skip=1, header=FALSE) names(fst)[c(3,5)] <- c("midp", "fst") plot(fst$midp, fst$fst, xlab="Physical position", ylab="Fst", col="#5f9ea0", pch=16) ``` --- # General feature format (GFF) from EnsemblPlants Maize [reference genome](https://plants.ensembl.org/Zea_mays/Info/Index) change to `largedata\lab4` folder: ```bash wget ftp://ftp.ensemblgenomes.org/pub/plants/release-46/fasta/zea_mays/dna/Zea_mays.B73_RefGen_v4.dna.chromosome.Mt.fa.gz ### then unzip it gunzip Zea_mays.B73_RefGen_v4.dna.chromosome.Mt.fa.gz ``` Similarly, we will download and unzip the [Mt GFF3](ftp://ftp.ensemblgenomes.org/pub/plants/release-46/gff3/zea_mays/Zea_mays.B73_RefGen_v4.46.chromosome.Mt.gff3.gz) file -- #### Use R to process the GFF3 file ```r # install.package("data.table") library("data.table") ## simply read in wouldn't work gff <- fread("largedata/lab4/Zea_mays.B73_RefGen_v4.46.chromosome.Mt.gff3", skip="#", header=FALSE, data.table=FALSE) ## grep -v means select lines that not matching any of the specified patterns gff <- fread(cmd='grep -v "#" largedata/lab4/Zea_mays.B73_RefGen_v4.46.chromosome.Mt.gff3', header=FALSE, data.table=FALSE) ``` --- # General feature format (GFF) version 3 ``` V1 V2 V3 V4 V5 V6 V7 V8 1 Mt Gramene chromosome 1 569630 . . . 2 Mt ensembl gene 6391 6738 . + . 3 Mt NCBI mRNA 6391 6738 . + . 4 Mt NCBI exon 6391 6738 . + . 5 Mt NCBI CDS 6391 6738 . + 0 6 Mt ensembl gene 6951 8285 . + . V9 1 ID=chromosome:Mt;Alias=AY506529.1,NC_007982.1;Is_circular=true 2 ID=gene:ZeamMp002;biotype=protein_coding;description=orf115-a1; 3 ID=transcript:ZeamMp002;Parent=gene:ZeamMp002; 4 Parent=transcript:ZeamMp002;Name=ZeamMp002.exon1;constitutive=1;ensembl_end_phase=0; 5 ID=CDS:ZeamMp002;Parent=transcript:ZeamMp002; 6 ID=gene:ZeamMp003;biotype=protein_coding;description=orf444 ``` --- # General feature format (GFF) version 3 ``` V1 V2 V3 V4 V5 V6 V7 V8 1 Mt Gramene chromosome 1 569630 . . . 2 Mt ensembl gene 6391 6738 . + . V9 1 ID=chromosome:Mt;Alias=AY506529.1,NC_007982.1;Is_circular=true 2 ID=gene:ZeamMp002;biotype=protein_coding;description=orf115-a1; ``` -------------- - 1 __sequence__: The name of the sequence where the feature is located. - 2 __source__: Keyword identifying the source of the feature, like a program (e.g. RepeatMasker) or an organization (like ensembl). - 3 __feature__: The feature type name, like "gene" or "exon". - In a well structured GFF file, all the children features always follow their parents in a single block (so all exons of a transcript are put after their parent "transcript" feature line and before any other parent transcript line). - 4 __start__: Genomic start of the feature, with a 1-base offset. - This is in contrast with other 0-offset half-open sequence formats, like [BED](). --- # General feature format (GFF) version 3 ``` V1 V2 V3 V4 V5 V6 V7 V8 1 Mt Gramene chromosome 1 569630 . . . 2 Mt ensembl gene 6391 6738 . + . V9 1 ID=chromosome:Mt;Alias=AY506529.1,NC_007982.1;Is_circular=true 2 ID=gene:ZeamMp002;biotype=protein_coding;description=orf115-a1; ``` -------------- - 5 __end__: Genomic end of the feature, with a 1-base offset. - This is the same end coordinate as it is in 0-offset half-open sequence formats, like BED. - 6 __score__: Numeric value that generally indicates the confidence of the source in the annotated feature. - A value of "." (a dot) is used to define a null value. - 7 __strand__: Single character that indicates the strand of the feature. - it can assume the values of "+" (positive, or 5' -> 3'), - "-", (negative, or 3' -> 5'), "." (undetermined). --- # General feature format (GFF) version 3 ``` V1 V2 V3 V4 V5 V6 V7 V8 1 Mt Gramene chromosome 1 569630 . . . 2 Mt ensembl gene 6391 6738 . + . V9 1 ID=chromosome:Mt;Alias=AY506529.1,NC_007982.1;Is_circular=true 2 ID=gene:ZeamMp002;biotype=protein_coding;description=orf115-a1; ``` -------------- - 8 __phase__: phase of CDS (__means CoDing Sequence__) features. - The phase indicates where the feature begins with reference to the reading frame. - The phase is one of the integers 0, 1, or 2, indicating the number of bases that should be removed from the beginning of this feature to reach the first base of the next codon. - 9 __attributes__: All the other information pertaining to this feature. - The format, structure and content of this field is the one which varies the most between the three competing file formats. --- # Work with GFF ```r names(gff) <- c("seq", "source", "feature", "start", "end", "score", "strand", "phase", "att") table(gff$feature) ``` ### Get genes and upstream and downstream 5kb regions ```r g <- subset(gff, feature %in% "gene") g$geneid <- gsub(".*gene:|;biotype.*", "", g$att) ### + strand gp <- subset(g, strand %in% "+") # nrow(gp) 75 ### get the 5k upstream of the + strand gene model gp_up <- gp gp_up$end <- gp_up$start - 1 gp_up$start <- gp_up$end - 5000 ### get the 5k downstream of the + strand gene model gp_down <- gp gp_down$start <- gp_down$end + 1 gp_down$end <- gp_down$start + 5000 ``` --- ### Get genes and upstream and downstream 5kb regions ```r ### - strand gm <- subset(g, strand %in% "-") dim(gm) # 82 fwrite(g, "cache/mt_gene.txt", sep="\t", row.names = FALSE, quote=FALSE) ``` --- ## Intepret the theta results ```r library("data.table") library("GenomicRanges") library("plyr") theta <- fread("cache/theta.txt", data.table=FALSE) names(theta)[1] <- "seq" up5k <- read.table("cache/mt_gene_up5k.txt", header=TRUE) ### define the subject file for theta values grc <- with(theta, GRanges(seqnames=seq, IRanges(start=Pos, end=Pos))) ### define the query file for genomic feature grf <- with(up5k, GRanges(seqnames=seq, IRanges(start=start, end=end), geneid=geneid)) ### find overlaps between the two tb <- findOverlaps(query=grf, subject=grc) tb <- as.matrix(tb) out1 <- as.data.frame(grf[tb[,1]]) out2 <- as.data.frame(grc[tb[,2]]) ### for each genomic feature, find the sites with non-missing data out <- cbind(out1, out2[, "start"]) names(out)[ncol(out)] <- "pos" ``` --- ## Intepret the theta results ```r #define unique identifier and merge with the thetas out$uid <- paste(out$seqnames, out$pos, sep="_") theta$uid <- paste(theta$seq, theta$Pos, sep="_") df <- merge(out, theta[, c(-1, -2)], by="uid") # for each upstream 5k region, how many theta values mx <- ddply(df, .(geneid), summarise, Pairwise = mean(Pairwise, na.rm=TRUE), thetaH = mean(thetaH, na.rm=TRUE), nsites = length(uid)) ``` --- ## Intepret the theta results Copy and paste everything above, and pack it into an `R` function: ```r get_mean_theta <- function(gf_file="cache/mt_gene_up5k.txt"){ # gf_file: gene feature file [chr, ="cache/mt_gene_up5k.txt"] theta <- fread("cache/theta.txt", data.table=FALSE) names(theta)[1] <- "seq" up5k <- read.table(gf_file, header=TRUE) ### define the subject file for theta values grc <- with(theta, GRanges(seqnames=seq, IRanges(start=Pos, end=Pos))) ### define the query file for genomic feature grf <- with(up5k, GRanges(seqnames=seq, IRanges(start=start, end=end), geneid=geneid)) ### find overlaps between the two tb <- findOverlaps(query=grf, subject=grc) tb <- as.matrix(tb) out1 <- as.data.frame(grf[tb[,1]]) out2 <- as.data.frame(grc[tb[,2]]) ### for each genomic feature, find the sites with non-missing data out <- cbind(out1, out2[, "start"]) names(out)[ncol(out)] <- "pos" #define unique identifier and merge with the thetas out$uid <- paste(out$seqnames, out$pos, sep="_") theta$uid <- paste(theta$seq, theta$Pos, sep="_") df <- merge(out, theta[, c(-1, -2)], by="uid") # for each upstream 5k region, how many theta values mx <- ddply(df, .(geneid), summarise, Pairwise = mean(Pairwise, na.rm=TRUE), thetaH = mean(thetaH, na.rm=TRUE), nsites = length(uid)) return(mx) } ``` --- ## Plot the results Run the customized R function ```r ### apply the function up5k <- get_mean_theta(gf_file="cache/mt_gene_up5k.txt") down5k <- get_mean_theta(gf_file="cache/mt_gene_down5k.txt") ``` And then plot the results: ```r library("ggplot2") up5k$feature <- "up 5k" down5k$feature <- "down 5k" res <- rbind(up5k, down5k) ggplot(res, aes(x=feature, y=Pairwise, fill=feature)) + geom_violin(trim=FALSE)+ labs(title="Theta value", x="", y = "Log10 (theta)")+ geom_boxplot(width=0.1, fill="white")+ scale_fill_brewer(palette="Blues") + theme_classic() ```