重测序分析脚本

重测序与群体遗传学
云南农业大学 云南生物资源保护与利用国家重点实验室
（李详 2019年12月31日）

软件安装

该部分软件很难用conda直接安装，安装步骤比较特殊。

lumpy-sv

1
2
3
4
5
6
7
lumpy-sv 依赖 Python2.7，切换到py27环境下进行安装。
conda activate py27
git clone --recursive https://github.com/arq5x/lumpy-sv.git
cd lumpy-sv
make
## 安装svtyper
conda install svtyper

cnvnator

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
cnvnator依赖于root软件包及samtools软件包(包含HTSlib)，因此先安装Samtools及root软件包。

## 安装samtools
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -xvf samtools-1.9.tar.bz2
cd samtools-1.9/
./configure
make
## 安装root软件, 直接下载后解压即可
wget https://root.cern/download/root_v6.18.04.Linux-ubuntu18-x86_64-gcc7.4.tar.gz
tar -zxvf root_v6.18.04.Linux-ubuntu18-x86_64-gcc7.4.tar.gz
## 安装cnvnator
git clone https://github.com/abyzovlab/CNVnator.git
cd CNVnator
## 链接samtools软件目录到cnvnator目录
ln -s ../samtools-1.9 ./samtools
## 加载root软件到环境变量
export ROOTSYS=/pub/software/root
export PATH=/pub/software/root/bin:$PATH
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib
## 安装cnvnator
make LIBS="-lcrypto" 

less ../data/sample.list | awk '{ print "bwa mem -t 2 -R '\''@RG\\tID:"$1"\\tSM:"$1"\\tPL:illumina'\'' ../01.ref/genome.fasta  ../data/"$1"_1.fq.gz ../data/"$1"_2.fq.gz |  /pub/software/samtools/samtools sort -@ 2 -m 1G -o "$1".sort.bam -" }'

变异检测

构建基因组index

1
2
3
4
5
6
7
8
9
10
11
# create link
ln -s  ../data/genome.fasta  ./genome.fasta

# for gatk
samtools faidx genome.fasta 

# for bwa
bwa index  genome.fasta  

# for picard
picard CreateSequenceDictionary R=genome.fasta

比对排序去重

1
2
3
4
5
6
7
8
9
10
11
# 比对排序
bwa mem -t 2 -R '@RG\tID:S1\tSM:S1\tPL:illumina' ../01.ref/genome.fasta  ../data/S1_1.fq.gz ../data/S1_2.fq.gz |  /pub/software/samtools/samtools sort -@ 2 -m 1G -o S1.sort.bam -
bwa mem -t 2 -R '@RG\tID:S2\tSM:S2\tPL:illumina' ../01.ref/genome.fasta  ../data/S2_1.fq.gz ../data/S2_2.fq.gz |  /pub/software/samtools/samtools sort -@ 2 -m 1G -o S2.sort.bam -

# 去除重复
picard -Xmx4g  MarkDuplicates I=S1.sort.bam O=S1.sort.rmdup.bam  CREATE_INDEX=true  REMOVE_DUPLICATES=true M=S1.marked_dup_metrics.txt
picard -Xmx4g  MarkDuplicates I=S2.sort.bam O=S2.sort.rmdup.bam  CREATE_INDEX=true  REMOVE_DUPLICATES=true M=S2.marked_dup_metrics.txt

# 比对率统计
/pub/software/samtools/samtools  flagstat  S1.sort.bam > S1.sort.bam.flagstat
/pub/software/samtools/samtools  flagstat  S2.sort.bam > S2.sort.bam.flagstat

SNP、Indel检测

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# step1.HaplotypeCaller
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" HaplotypeCaller -R ../01.ref/genome.fasta -I ../02.mapping/S1.sort.rmdup.bam  -ERC GVCF -O S1.g.vcf 1>S1.HC.log   2>&1
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" HaplotypeCaller -R ../01.ref/genome.fasta -I ../02.mapping/S2.sort.rmdup.bam  -ERC GVCF -O S2.g.vcf 1>S2.HC.log   2>&1

# step2.CombineGVCFs
ls ./S*.g.vcf > gvcf.list
gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"    CombineGVCFs -R ../01.ref/genome.fasta -V gvcf.list  -O all.merge.g.vcf 

# step3.GenotypeGVCFs
gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"   GenotypeGVCFs -R ../01.ref/genome.fasta --variant all.merge.g.vcf -O all.merge_raw.vcf

# step4.SNP
gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  SelectVariants  -R ../01.ref/genome.fasta -V all.merge_raw.vcf --select-type SNP -O all.raw.snp.vcf

gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  VariantFiltration -R ../01.ref/genome.fasta -V all.raw.snp.vcf --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name 'SNP_filter' -O all.filter.snp.vcf

gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  SelectVariants  -R ../01.ref/genome.fasta -V all.filter.snp.vcf --exclude-filtered  -O all.filtered.snp.vcf

# step5.InDel
gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  SelectVariants  -R ../01.ref/genome.fasta -V all.merge_raw.vcf --select-type INDEL -O all.raw.indel.vcf

gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  VariantFiltration -R ../01.ref/genome.fasta -V all.raw.indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name 'INDEL_filter' -O all.filter.indel.vcf

gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  SelectVariants  -R ../01.ref/genome.fasta -V all.filter.indel.vcf  --exclude-filtered   -O all.filtered.indel.vcf

SV检测

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# step1.discordants_bam
/pub/software/samtools/samtools view -b -F 1294 -@ 12 ../02.mapping/S1.sort.rmdup.bam > S1.discordants.bam
/pub/software/samtools/samtools view -b -F 1294 -@ 12 ../02.mapping/S2.sort.rmdup.bam > S2.discordants.bam

# step2.splitters_bam
/pub/software/samtools/samtools view -h ../02.mapping/S1.sort.rmdup.bam | /pub/software/lumpy-sv/scripts/extractSplitReads_BwaMem -i stdin | /pub/software/samtools/samtools  sort -   > S1.splitters.bam
/pub/software/samtools/samtools view -h ../02.mapping/S2.sort.rmdup.bam | /pub/software/lumpy-sv/scripts/extractSplitReads_BwaMem -i stdin | /pub/software/samtools/samtools  sort -   > S2.splitters.bam

# step3.lumpy
/pub/software/lumpy-sv/bin/lumpyexpress -B ../02.mapping/S1.sort.rmdup.bam,../02.mapping/S2.sort.rmdup.bam  -S S1.splitters.bam,S2.splitters.bam   -D S1.discordants.bam,S2.discordants.bam -o all.sv.lumpy.vcf

# step4.ind_genotype
vcftools --vcf all.sv.lumpy.vcf --indv S1 --recode --recode-INFO-all --out S1 && svtyper -i S1.recode.vcf -B ../02.mapping/S1.sort.rmdup.bam  -o S1.genotype.vcf
vcftools --vcf all.sv.lumpy.vcf --indv S2 --recode --recode-INFO-all --out S2 && svtyper -i S2.recode.vcf -B ../02.mapping/S2.sort.rmdup.bam  -o S2.genotype.vcf

CNV检测

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# step1.prepare_genome

## split genome into single sequence

 ln -s  ../data/genome.fasta ./genome.fa

 seqkit split -i ./genome.fa  -O split

 ls ./split/genome.id_*.fa |sed 's/\.\/split\/genome\.id_//' |awk '{print "mv ./split/genome.id_"$aa"  ./split/"$aa}' > mv_name.sh

 sh mv_name.sh

# step2.cnvnator
export ROOTSYS=/pub/software/root
export PATH=/pub/software/root/bin:$PATH
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib

## 链接文件到当前目录
ln -s ../02.mapping/S1.sort.rmdup.bam 

## 提取比对结果
/pub/software/CNVnator/cnvnator -genome genome.fa -root S1.root -tree S1.sort.rmdup.bam

##生成深度分布
/pub/software/CNVnator/cnvnator -genome genome.fa -root S1.root -his 500  -d split

## 进行统计计算
/pub/software/CNVnator/cnvnator -root S1.root -stat 500

## 检查bin size是否合适
/pub/software/CNVnator/cnvnator -root S1.root -eval 500  > S1.eval.ratio 

## RD信号分割
/pub/software/CNVnator/cnvnator -root S1.root -partition 500

## 进行CNV检测
/pub/software/CNVnator/cnvnator -root S1.root  -call 500 > S1.cnv

SNP INDEL注释

1
2
3
4
5
6
7
8
9
10
11
12
13
# step1.prepare_genome
ln -s ../data/genome.gtf 
ln -s ../data/genome.fasta

gtfToGenePred  -genePredExt genome.gtf genome_refGene.txt
/pub/software/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile genome.fasta genome_refGene.txt --out genome_refGeneMrna.fa

#step2.ann_vcf
ln -s  ../03.SNP_indel/all.filtered.snp.vcf 

/pub/software/annovar/convert2annovar.pl -format vcf4old all.filtered.snp.vcf  >all.snp.vcf.annovar.input

/pub/software/annovar/annotate_variation.pl -geneanno --neargene 2000 -buildver genome -dbtype refGene -outfile all.anno -exonsort all.snp.vcf.annovar.input ./

CNV和SV注释

1
2
3
4
5
6
7
less -S ../04.SV/all.sv.lumpy.vcf|awk -F ";" '{ if($1!~/#/){print $1"\t"$2"\t"$3"\t"$4} }'|awk '{print $1"\t"$2"\t"$11"\t"$3"\t"$8}'|sed 's/END=//' |sed 's/SVTYPE=//'| awk '$3<$2' > SV.bed
less -S ../04.SV/all.sv.lumpy.vcf|awk -F ";" '{ if($1!~/#/){print $1"\t"$2"\t"$3"\t"$4} }'|awk '{print $1"\t"$2"\t"$11"\t"$3"\t"$8}'|sed 's/END=//' |sed 's/SVTYPE=//'  > SV.bed
awk '$5 !~ /BND/' SV.bed >  SV_filter.bed
bedtools  intersect -wo  -a SV_filter.bed -b gene.gtf  > SV_filter.Anno
less -S ../05.CNV/S1.cnv | awk '{print $2 "\t" $1}'|sed 's/:/\t/'|sed 's/-/\t/' > CNV.bed
bedtools  intersect -wo  -a  CNV.bed  -b gene.gtf > CNV.Anno

群体结构

过滤

1
2
3
4
5
6
7
8
9
vcf=../data/all.vcf

## filter missing and maf
plink --vcf  $vcf --geno 0.1 --maf 0.01 --out all.missing_maf --recode vcf-iid  --allow-extra-chr --set-missing-var-ids @:# --keep-allele-order

## filter LD 
plink --vcf  all.missing_maf.vcf  --indep-pairwise 50 10 0.2 --out tmp.ld   --allow-extra-chr --set-missing-var-ids @:# 
plink --vcf  all.missing_maf.vcf  --make-bed --extract tmp.ld.prune.in  --out all.LDfilter --recode vcf-iid  --keep-allele-order  --allow-extra-chr --set-missing-var-ids @:#

建树

NJ

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
perl  ../../script/vcf2phy.pl  ../../00.filter/all.LDfilter.vcf  > sequences.dat

echo -e "sequences.dat\nY" > dnadist.cfg

dnadist < dnadist.cfg  >dnadist.log

mv outfile infile.dist

echo -e "infile.dist\nY"  > neighbor.cfg

neighbor  <  neighbor.cfg  >nj.log

less infile.dist | tr '\n' '|'| sed 's/| / /g' | tr '|' '\n' >infile.dist.table
less outtree | tr '\n' ' '|sed 's/ //g' > outtree.nwk

tree_snphylo

1
2
/pub/software/SNPhylo/snphylo.sh -v ../../00.filter/all.LDfilter.vcf -r -l 1 -m 0.01 -o GA0001 -b -B 2 -P tree_snphylo

PCA

1
2
3
4
plink --vcf  ../00.filter/all.LDfilter.vcf --pca 10 --out  PCA_out   --allow-extra-chr --set-missing-var-ids @:# 

Rscript ../script/draw_PCA.R  PCA_out.eigenvec 1 2  ../data/sample.pop  PCA_out_figure

群体结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# step0.vcf2bed
plink --vcf  ../00.filter/all.LDfilter.vcf  --make-bed --out all  --allow-extra-chr --keep-allele-order --set-missing-var-ids @:#

# step1.get_Admix_shell
seq 2 4 | awk '{print "admixture --cv -j2 all.bed "$1" 1>admix."$1".log 2>&1"}' > admixture.sh

# step2.run_admix
sh admixture.sh

# step3.draw_structure
mkdir result
cp  ./*.Q result/
Rscript  ../script/draw_admixture.R result all.nosex  structure 

# admixture
admixture --cv -j2 all.bed 2 1>admix.2.log 2>&1
admixture --cv -j2 all.bed 3 1>admix.3.log 2>&1
admixture --cv -j2 all.bed 4 1>admix.4.log 2>&1

遗传衰减

1
2
3
4
5
6
7
8
9
10
# step1.LD_calcu
/pub/software/PopLDdecay/bin/PopLDdecay -InVCF  ../data/all.vcf -SubPop  ../data/pop.SC.table -MaxDist 500 -OutStat pop.SC.stat
/pub/software/PopLDdecay/bin/PopLDdecay -InVCF  ../data/all.vcf -SubPop ../data/pop.YZR.table -MaxDist 500 -OutStat pop.YZR.stat

# step2.draw_singlePlot
perl /pub/software/PopLDdecay/bin/Plot_OnePop.pl -inFile pop.SC.stat.gz -output pop.SC.ld
perl /pub/software/PopLDdecay/bin/Plot_OnePop.pl -inFile pop.YZR.stat.gz -output pop.YZR.ld

# step3.draw_multiPlotare/PopLDdecay/bin/Plot_MultiPop.pl -inList ld_stat.list -output ld_stat.multi

选择分析

i_ROD

1
2
3
4
5
6
7
8
9
10
11
12
# step1.pi
vcf=../../data/all.vcf
window=20000
step=2000
vcftools  --vcf $vcf --window-pi $window --window-pi-step  $step  --keep ../../data/pop.SC.table   --out ./Pi.pop1
vcftools  --vcf $vcf --window-pi $window --window-pi-step  $step  --keep ../../data/pop.YZR.table  --out ./Pi.pop2

# step2.ROD
pi1=Pi.pop1.windowed.pi
pi2=Pi.pop2.windowed.pi
perl ../../script/merge2pi_ROD.pl $pi1 $pi2  >Pi.ROD

tajimaD

1
2
3
4
5
6
vcf=../../data/all.vcf
window=20000

vcftools --vcf $vcf --TajimaD  $window  --keep  ../../data/pop.SC.table  --out TajimaD.pop1
vcftools --vcf $vcf --TajimaD  $window  --keep  ../../data/pop.YZR.table --out   TajimaD.pop2

Fst

1
2
3
4
5
vcf=../../data/all.vcf
window=20000
step=2000
vcftools  --vcf $vcf --fst-window-size $window --fst-window-step $step  --weir-fst-pop  ../../data/pop.SC.table --weir-fst-pop ../../data/pop.YZR.table --out  ./Fst.pop1.pop2

ROD_Fst

1
2
3
4
5
6
7
8
9
10
11
12
13
fst=../03.Fst/Fst.pop1.pop2.windowed.weir.fst
ROD=../01.pi_ROD/Pi.ROD
gff=../../data/genome.gff

perl ../../script/top_Fst_ROD_S.pl  $fst $ROD 0.05 0.05  Fst_ROD.table Fst_ROD.stat

awk '$NF=="top"' Fst_ROD.table  > Fst_ROD.table.top


awk '$3=="gene"' $gff > gene.gff

bedtools  intersect -wo -F 0.1  -a  Fst_ROD.table.top -b gene.gff | awk -F "\t"  '{print $7"\t"$10"\t"$11"\t"$13"\t"$15}' | sort -u  > Fst_ROD.table.top.gene

GWAS

mpute

1
2
beagle -Xmx4g  -Djava.io.tmpdir=./TMP  gt=../data/all.vcf   out=all.impute impute=true window=10 nthreads=2

GWAS

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# step1.prepare
## prepare vcf
plink --vcf ../data/all.vcf --maf 0.05 --geno 0.1  --recode12  --output-missing-genotype 0 --transpose --out snp_filter   --set-missing-var-ids @:#  --allow-extra-chr


## prepare phenotype

perl ../script/sort_pheno.pl snp_filter.tfam ../data/trait.table > trait.sort.txt

# step2.kinship
/pub/software/emmax/emmax-kin-intel64 -v -d 10  -o ./pop.kinship  snp_filter

# step3.gwas
/pub/software/emmax/emmax-intel64 -v -d 10 -t snp_filter  -p trait.sort.txt -k pop.kinship   -o emmax.out 1> emmax.log 2>emmax.err

# draw_manhattan
paste  snp_filter.map  emmax.out.ps | awk  'BEGIN{print "SNP\tCHR\tBP\tP"}{if($2==$5){print $2"\t"$1"\t"$4"\t"$NF}}'  > emmax.out.ps.manht_input
Rscript ../script/manhattan.R emmax.out.ps.manht_input  emmax.out.ps.manht_figure

GWAS-Q

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# step3.gwas
ln -s ../01.GWAS/snp_filter.tped
ln -s ../01.GWAS/snp_filter.tfam
ln -s ../01.GWAS/snp_filter.nosex
ln -s ../01.GWAS/snp_filter.map

ln -s ../01.GWAS/trait.sort.txt
ln -s ../../02.population_genetics/03.structure/all.3.Q

paste snp_filter.nosex all.3.Q |awk '{print $1" "$1" 1 "$3" "$4}'  > pop.Qmatrix

/pub/software/emmax/emmax-intel64 -v -d 10 -t snp_filter  -p trait.sort.txt -k pop.kinship -c pop.Qmatrix   -o emmax.out 1> emmax.log 2>emmax.err

# step4.draw_manhattan
paste  snp_filter.map  emmax.out.ps | awk  'BEGIN{print "SNP\tCHR\tBP\tP"}{if($2==$5){print $2"\t"$1"\t"$4"\t"$NF}}'  > emmax.out.ps.manht_input
Rscript ../script/manhattan.R emmax.out.ps.manht_input  emmax.out.ps.manht_figure_Q

GWAS_realData

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# step1.prepare
## prepare vcf
plink --vcf  ../data_real/Sample215.M5M8H3.2allel.vcf  --maf 0.05 --geno 0.1  --recode12  --output-missing-genotype 0 --transpose --out snp_filter   --set-missing-var-ids @:#  --allow-extra-chr


## prepare phenotype

perl ../script/sort_pheno.pl snp_filter.tfam  ../data_real/trait.C16_0.table  > trait.sort.txt

# step2.kinship
/pub/software/emmax/emmax-kin-intel64 -v -d 10  -o ./pop.kinship  snp_filter

# step3.gwas
/pub/software/emmax/emmax-intel64 -v -d 10 -t snp_filter  -p trait.sort.txt -k pop.kinship   -o emmax.out 1> emmax.log 2>emmax.err

# step4.draw_manhattan
paste  snp_filter.map  emmax.out.ps | awk  'BEGIN{print "SNP\tCHR\tBP\tP"}{if($2==$5){print $2"\t"$1"\t"$4"\t"$NF}}'  > emmax.out.ps.manht_input
Rscript ../script/manhattan.R emmax.out.ps.manht_input  emmax.out.ps.manht_figure_readData

QTL-seq

change_VCF2Table

1
2
3
4
5
ref=
vcf=
outfile=
gatk  VariantsToTable -R  $ref -V $vcf -F CHROM -F POS -F REF -F ALT -GF AD -GF DP -GF GQ -GF PL -O  $outfile

QTLSeq_Analysis

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
library(QTLseqr)
library("ggplot2")


HighBulk <- "SRR834931" ## 高值混池名称
LowBulk <- "SRR834927"  ## 低值混池名称
Chroms <- paste0(rep("Chr", 12), 1:12) ## 染色体列表

## SNP数据读取
df <- importFromGATK(file = "../data/Yang_et_al_2013.table" ,
    highBulk = HighBulk, # 指定高值混池
    lowBulk = LowBulk, # 指定低至混池
    chromList = Chroms # 指定分析用染色体列表
    ) 


## 绘制深度分布图 

p1 <- ggplot(data = df) +
  geom_histogram(aes(x = DP.HIGH + DP.LOW)) +
  xlim(0,800)
  
pdf(file="SNP_depth.pdf")
p1
dev.off()

## 绘制ref等位基因频率分布图
p2 <- ggplot(data = df) +
  geom_histogram(aes(x = REF_FRQ))
  
pdf(file="ref_allele_frequency.pdf")
p2
dev.off()


## 绘制高值混池SNP-index分布
p3 <- ggplot(data = df) +
  geom_histogram(aes(x = SNPindex.HIGH))
pdf(file="SNPindex.HIGH.dis.pdf")
p3
dev.off()

## 绘制低值混池SNP-index分布
p4 <- ggplot(data = df) +
  geom_histogram(aes(x = SNPindex.LOW))

pdf(file="SNPindex.LOW.dis.pdf")
p4
dev.off()

## SNP过滤

df_filt <- filterSNPs(
    SNPset = df,
    refAlleleFreq = 0.20, ## ref allele频率过滤，0.2表示0.2~0.8之间
    minTotalDepth = 100, ## 最小深度过滤
    maxTotalDepth = 400, ## 最大深度过滤
    minSampleDepth = 40, ## 单个样品最小深度
    minGQ = 99, ## genotype quality 过滤
    verbose = TRUE  ## 输出日志
  )

## 进行deltaSNPindex 计算
df_filt <- runQTLseqAnalysis(df_filt,
    windowSize = 1e6, ## 窗口大小
    popStruc = "F2",  ## 群体类型，F2 或者RIL
    bulkSize = c(385, 430), ## 混样个数，第一个高值组
    replications = 10000, ## bootstrap次数
    intervals = c(95, 99) ## 置信区间
)


## SNP 沿染色体分布图
p5 <- plotQTLStats(SNPset = df_filt, var = "nSNPs")
pdf(file="SNP_filter.window.pdf", width = 20, height = 4)
p5
dev.off()

## deltaSNPindex沿染色体分布
p6 <- plotQTLStats(SNPset = df_filt, var = "deltaSNP", plotIntervals = TRUE)
pdf(file="deltaSNPindex.pdf", width = 20, height = 4)
p6
dev.off()

## 提取显著性区域
QTL <- getSigRegions(SNPset = df_filt, 
    method = "QTLseq", interval = 95)

## 输出到文件
write.table(QTL[[1]],
    "QTLseq_result.SigRegions.table",
    sep="\t",quote=F)

## 提取QTLseq结果
results99 <- getQTLTable(
  SNPset = df_filt,  
  method = "QTLseq", 
  interval = 99,  ## 结果阈值
  export = FALSE)

## 输出到文件
write.table(results99,
    file="QTLseq_result_deltaSNPindex_CI99.table" ,
    sep="\t", quote = F)

💌lixiang117423@gmail.com

💌lixiang117423@foxmail.com

生物信息学 - 群体遗传学 - Linux
#生物信息学 #Linux
重测序分析脚本
https://lixiang117423.github.io/article/5269tghy/
作者
小蓝哥
发布于
2021年11月23日
许可协议