目錄
- 前言
- 方法簡介
- 從vcf文件里提取有用信息
- tidy vcfR
- vcf可視化1
- vcf可視化2
- 測序深度覆蓋度
- 窗口縮放
- 如何單獨(dú)分離染色體
- 利用vcf信息判斷物種染色體倍數(shù)
- CNV分析
1.讀取數(shù)據(jù)
vcfR可以直接讀取vcf格式的數(shù)據(jù)。如果同時(shí)讀取參照序列fasta格式的序列文件和gff格式文件的注釋文件還可以獲取更完整的信息(此步驟并非必須，可以只讀取vcf數(shù)據(jù))带斑。在此處便于重復(fù)用到了pinfsc50包碗殷。這個(gè)包里是植物致病微生物的基因序列測序結(jié)果孤紧。包含了一個(gè)vcf文件誊册，一個(gè)fasta文件和一個(gè)gff文件仙蛉。
1.1 整合數(shù)據(jù)
pkg <- "pinfsc50"
vcf_file <- system.file("extdata", "pinf_sc50.vcf.gz", package = pkg)
dna_file <- system.file("extdata", "pinf_sc50.fasta", package = pkg)
gff_file <- system.file("extdata", "pinf_sc50.gff", package = pkg)
1.2 用vcfR讀取vcf數(shù)據(jù)
library(vcfR)
vcf <- read.vcfR( vcf_file, verbose = FALSE )
1.3 用ape讀取fasta數(shù)據(jù)
這里用到參照序列的數(shù)據(jù)材彪。
dna <- ape::read.dna(dna_file, format = "fasta")
1.4 讀取gff格式的注釋文件
gff <- read.table(gff_file, sep="\t", quote="")
當(dāng)這些數(shù)據(jù)被讀取到內(nèi)存的時(shí)候就可以開始對染色體名字或者其它一些東西進(jìn)行修改了。由于vcfR更擅長對的單獨(dú)染色體進(jìn)行分析觉痛，所以當(dāng)你的基因過大或者有很多樣本的時(shí)候役衡，建議對數(shù)據(jù)進(jìn)行拆分。
2. 建立chromR項(xiàng)目
讀取完數(shù)據(jù)以后就可以建立chromR薪棒，來對數(shù)據(jù)進(jìn)行詳細(xì)的分析手蝎。
library(vcfR)
chrom <- create.chromR(name='Supercontig', vcf=vcf, seq=dna, ann=gff)
3. 對chromR進(jìn)行進(jìn)一步處理
首先對數(shù)據(jù)進(jìn)行初步的可視化,
plot(chrom)
我們在上面的圖里得到很多信息，比方說測序深度(DP)的峰在500俐芯，但是拖著尾巴棵介，這個(gè)尾巴表示數(shù)據(jù)里包含著CNV信息。然后比對質(zhì)量(MQ)的峰值在60,于是我們可以以60為中心對數(shù)據(jù)進(jìn)行過濾吧史。
使用masker可以對數(shù)據(jù)進(jìn)行過濾標(biāo)記邮辽。然后可視化過濾以后的數(shù)據(jù)。
chrom <- masker(chrom, min_QUAL = 1, min_DP = 300, max_DP = 700, min_MQ = 59.9, max_MQ = 60.1)
plot(chrom)
是不是順眼多了贸营。當(dāng)然我們也可以看一下SNP的分布情況吨述。注意右下角的圖。
chrom <- proc.chromR(chrom, verbose=TRUE)
plot(chrom)
4. 數(shù)據(jù)可視化
用chromoqc()可以對數(shù)據(jù)進(jìn)行更完整的可視化钞脂。包括外顯子內(nèi)含子的分布揣云，GC含量的分布等等。
chromoqc(chrom, dp.alpha=20)
5. 數(shù)據(jù)輸出
最后可以用函數(shù)write.vcf()把數(shù)據(jù)輸出成新的vcf文件芳肌。
之后的文章里會(huì)有更加詳細(xì)的說明
