目录
从参照序列提取单独染色体
会用到ape包，这里以提取染色体"Supercontig_1.1"为例。
library(ape)
dna <- read.dna("pinf_super_contigs.fa", format = "fasta")
dna2 <- dna[ grep( "Supercontig_1.1 ", names(dna) ) ]
names(dna2) <- "Supercontig_1.1"
dna2 <- as.matrix(dna2)
注意一下，此处用了names()函数来重命名染色体。
从注释文件提取染色体
提取注释文件的信息其实也很简单。
gff <- read.table('gff_file.gff', sep="\t", quote="")
gff2 <- gff[grep("Supercontig_1.1", gff[,1]),]
从vcf文件提取染色体
注意，不是在R里完成的。需要在Terminal里进行。windows用户可以用vcftools。
grep "^#" my_variants.vcf > header.vcf # Meta
grep "^Supercontig_1.1" my_variants.vcf > tmp.vcf # Body
cat header.vcf tmp.vcf > sc1.vcf.gz
如果使用的是压缩文件的话，
zgrep "^#" my_variants.vcf.gz | gzip -c > header.vcf.gz # Meta
zgrep "^Supercontig_1.1" my_variants.vcf.gz > tmp.vcf.gz # Body
zcat header.vcf.gz tmp.vcf.gz | gzip -c > sc1.vcf.gz
可以解决问题。
使用vcfR
然后按照之前介绍的方法进行分析，可视化。
library(vcfR)
vcf <- read.vcfR("sc1.vcf")
chrom <- create.chromR(name='Supercontig', vcf=vcf, seq=dna2, ann=gff2)
