5. Write new VCFs - DianaCarolinaVergara/SNPs_pipeline GitHub Wiki

Write new filtered VCFs

After all the filtering processes we have to create a new vcf file with the new changes.

Example:

write.vcf(example_DP_filtered,file= "example_novo_DP_filtered.vcf")

example_novo_DP_filtered_vcf<- read.vcfR("example_novo_DP_filtered.vcf" , verbose = TRUE )
example_novo_DP_filtered_vcf
head(example_novo_DP_filtered_vcf)

And with Muricea, doing a violin plot to visualize the changes in the new vcf.

#Write new filtered VCFs

write.vcf(Muricea_DP_filtered,file= "Muricea_DP_90Miss_filtered.vcf")

violin plot code __

Muricea_DP_90Miss_filtered.vcf <- read.vcfR("Muricea_DP_90Miss_filtered.vcf")
Muricea_DP_90Miss_filtered.vcf

dp <- extract.gt(Muricea_DP_90Miss_filtered.vcf,  element = "DP", as.numeric = TRUE)
class(dp)


dpf <- melt(dp, varnames = c("Index", "Sample"),
            value.name = "Depth", na.rm = TRUE)
dpf <- dpf[ dpf$Depth > 0, ]
p <- ggplot(dpf, aes(x = Sample, y = Depth))
p <- p + geom_violin(fill = "#C0C0C0", adjust = 1.0,
                     scale = "count", trim = TRUE)
p <- p + theme_bw()
p <- p + theme(axis.title.x = element_blank(),
               axis.text.x = element_text(angle = 60, hjust = 1))
p <- p + scale_y_continuous(trans = scales::log2_trans(),
                            breaks = c(1, 10, 100, 800),
                            minor_breaks = c(1:10, 2:10 * 10, 2:8 * 100))
p <- p + theme(panel.grid.major.y = element_line(color = "#A9A9A9", size = 0.6))
p <- p + theme(panel.grid.minor.y = element_line(color = "#C0C0C0", size = 0.2))
p <- p + ylab("Depth (DP)")
p