简化基因组数据分析实战

mkdir -p stacks-exercisecd stacks-exercise
wget -q http://catchenlab.life.illinois.edu/data/rochette2017_gac_or.tar.gz
tar xf http://catchenlab.life.illinois.edu/data/rochette2017_gac_or.tar.gz

# 适用conda安装软件conda install -c bioconda bwa
conda install -c bioconda samtools
conda install -c bioconda stacks=1.47

mkdir -p stacks-exercise/{00-raw-data,01-clean-data,02-read-alignment,reference/genome,stacks/{de-novo,ref},test/{de-novo,ref},info}# 目录结构stacks-exercise/
|-- 00-raw-data # 原始数据|-- 01-clean-data # 处理后数据|-- 02-read-alignment # 比对后数据|-- info # barcode信息|-- reference # 参考基因组|   |-- genome
|-- rochette2017_gac_or.tar.gz
|-- stacks # 实际运行|   |-- de-novo
|   |-- ref
|-- test # 测试数据集|-- de-novo
|-- ref

cd stacks-exercise
mv rochette2017_gac_or/top/raw/{lane1,lane2,lane3} 00-raw-data

mv rochette2017_gac_or/top/info/*.tsv info/# barcode和样本的对应关系head -n3 info/barcodes.lane1.tsv
CTCGCC  sj_1819.35GACTCT  sj_1819.31GAGAGA  sj_1819.32# 样本和群体的对应关系head -n3 info/popmap.tsv
cs_1335.01  cs
cs_1335.02  cs
cs_1335.03  cs

<6-nt barcode>TGCAGG

# 在项目根目录下raw_dir=00-raw-data/lane1
barcodes_file=info/barcodes.lane1.tsv
process_radtags -p $raw_dir -b $barcode_file \
-o 01-clean-data/ -e sbfI --inline_null \
-c -q -r &> 01-clean-data/process_radtags.lane1.oe &

# 位于项目根目录下
cd reference/genome
wget -q ftp://ftp.ensembl.org/pub/release-91/fasta/gasterosteus_aculeatus/dna/Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa.gz
gzip -d Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa.gzcd ..
mkdir -p index/bwa/
genome_fa=genome/Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa
bwa index -p index/bwa/gac $genome_fa &> index/bwa/bwa_index.oe
# 结果如下
|-- genome
|   |-- Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa
|-- index
    |-- bwa
        |-- bwa_index.oe
        |-- gac.amb
        |-- gac.ann
        |-- gac.bwt
        |-- gac.pac
        |-- gac.sa

# 位于项目根目录下# 在测试文件下按照参数创建文件夹
mkdir -p test/ref/{alignment-bwa,stacks-bwa}
## bwa-mem比对
sample=cs_1335.01fq_file=01-clean-data/$sample.fq.gz
bam_file=test/ref/alignment-bwa/${sample}_default.bam
bwa_index=reference/index/bwa/gac
bwa mem -M $bwa_index $fq_file | samtools view -b > $bam_file &

samtools flagstat test/ref/alignment-bwa-default/cs_1335.01_default.bam
#1310139   0 in total (QC-passed reads   QC-failed reads)
#7972   0 secondary
#0   0 supplementary
#0   0 duplicates#1271894   0 mapped (97.08% : N/A)

# 位于项目根目录下
sample=cs_1335.01fq_file=01-clean-data/$sample.fq.gz
bam_file=test/ref/alignment-bwa/${sample}_B3_O55.bam
bwa_index=reference/index/bwa/gac
bwa mem -M -B 3 -O 5,5 $bwa_index $fq_file | samtools view -b > $bam_file &
samtools flagstat test/ref/alignment-bwa-default/cs_1335.01_default.bam
#1309830   0 in total (QC-passed reads   QC-failed reads)
#7663   0 secondary
#0   0 supplementary
#0   0 duplicates#1272297   0 mapped (97.13% : N/A)

# 位于项目根目录下
sample=cs_1335.01
sample_index=1
bam_file=test/ref/alignment-bwa/${sample}_B3_O55.bam
log_file=test/ref/stacks-bwa/$sample.pstacks.oe
pstacks -t bam -f $bam_file -i $sample_index -o test/ref/stacks-bwa/ &> $log_file

cstacks -p 10 --aligned -P 03-alignment-stacks/ -M info/popmap.tsv
# 以其中一个样本为例
sample=cs_1335.01
log_file=$sample.sstacks.oe
sstacks --aligned -c 03-alignment-stacks/batch_1 -s 03-alignment-stacks/$sample -o 03-alignment-stacks/ &> 03-alignment-stacks/$log_file

INDEX = "reference/index/bwa/gac"
 
SAMPLES, = glob_wildcards("01-raw-data/{sample}.fq.gz")
 
INDEX_DICT = {value: key for key, value in dict(enumerate(SAMPLES, start=1))}
 
FQ_FILES = expand("01-clean-data/{sample}.fq.gz", sample=SAMPLES)
 
BAM_FILES = expand("02-read-alignment/ref-based/{sample}.bam", sample=SAMPLES)
 
SNP_TSVS = expand("03-alignment-stacks/{sample}.snp.tsv.gz", sample=SAMPLES)
 
CATALOG = "03-alignment-stacks/batch_1.catalog.snps.tsv.gz"
 
MATCHES = expand("03-alignment-stacks/{sample}.matches.tsv.gz", sample=SAMPLES)
 
rule all:
 
 input:
 
 BAM_FILES,
 
 
rule bwa_mem:
 
 input: "01-clean-data/{sample}.fq.gz"
 
 params:
 
 index = INDEX,
 
 mismatch = "3",
 
 gap = "5,5"
 
 threads: 8
 
 output: "02-read-alignment/ref-based/{sample}.bam"
 
 shell:"""
 
 mkdir -p 02-read-alignment
 
 bwa mem -t {threads} -M -B {params.mismatch} -O {params.gap} {params.index} {input} \
 
 | samtools view -b > {output}
 
 """
 
 
 rule pstacks:
 
 input: "02-read-alignment/ref-based/{sample}.bam"
 
 params:
 
 index = lambda wildcards: INDEX_DICT.get(wildcards.sample),
 
 outdir = "03-alignment-stacks/"
 
 threads: 8
 
 output:
 
 "03-alignment-stacks/{sample}.snps.tsv.gz",
 
 "03-alignment-stacks/{sample}.tags.tsv.gz",
 
 "03-alignment-stacks/{sample}.models.tsv.gz",
 
 "03-alignment-stacks/{sample}.alleles.tsv.gz"
 
 log: "03-alignment-stacks/{sample}.pstacks.oe"
 
 shell:"""
 
 mkdir -p 03-alignment-stacks
 
 pstacks -p {threads} -t bam -f {input} -i {params.index} -o {params.outdir} &> {log}
 
 """
 
# A catalog can be built from any set of samples processed by the ustacks or pstacks programs
 
rule cstacks:
 
 input: "info/popmap.tsv",expand("03-alignment-stacks/{sample}.snps.tsv.gz", sample=SAMPLES)
 
 threads: 10
 
 output:
 
 "03-alignment-stacks/batch_1.catalog.alleles.tsv.gz",
 
 "03-alignment-stacks/batch_1.catalog.snps.tsv.gz",
 
 "03-alignment-stacks/batch_1.catalog.tags.tsv.gz"
 
 shell:
 
 "cstacks -p {threads} --aligned -P 03-alignment-stacks/ -M {input[0]}"
 
 
# Sets of stacks, i.e. putative loci, constructed by the ustacks or pstacks programs
 
# can be searched against a catalog produced by cstacks.
 
rule sstacks:
 
 input:
 
 "03-alignment-stacks/{sample}.snps.tsv.gz",
 
 "03-alignment-stacks/{sample}.tags.tsv.gz",
 
 "03-alignment-stacks/{sample}.models.tsv.gz",
 
 "03-alignment-stacks/{sample}.alleles.tsv.gz"
 
 params:
 
 catalog = "03-alignment-stacks/batch_1",
 
 sample = lambda wildcards: "03-alignment-stacks/"   wildcards.sample
 
 output:
 
 "03-alignment-stacks/{sample}.matches.tsv.gz"
 
 log: "03-alignment-stacks/{sample}.sstacks.oe"
 
 shell:"""
 
 sstacks --aligned -c {params.catalog} -s {params.sample} -o 03-alignment-stacks &> {log}
 
 """

snakemake --cluster "qsub -V -cwd" -j 20 --local-cores 10 &

cat popmap_sample.tsv
 
cs_1335.01 cs
cs_1335.02 cs
cs_1335.19 cs
pcr_1193.00 pcr
pcr_1193.01 pcr
pcr_1193.02 pcr
pcr_1213.02 pcr
sj_1483.01 sj
sj_1483.02 sj
sj_1483.03 sj
wc_1218.04 wc
wc_1218.05 wc
wc_1218.06 wc
wc_1219.01 wc

mkdir -p test/de-novo/stacks_M{1..9}

# 项目根目录
 
M=1
 
popmap=info/popmap_sample.tsv
 
reads_dir=01-clean-data
 
out_dir=test/de-novo/stacks_M${M}
 
log_file=$out_dir/denovo_map.oe
 
threads=8
 
denovo_map.pl -T ${threads} --samples ${reads_dir} -O ${popmap} -o ${out_dir} -M $M -n $M -m 3 -b 1 -S &> ${log_file} &

grep -iE "(err|e:|warn|w:|fail|abort)" test/de-novo/stacks_M1/denovo_map.oe

# 项目根目录
 
for M in `seq 1 9`
 
do
 
popmap=info/popmap_sample.tsv
 
stacks_dir=test/de-novo/stacks_M${M}
 
out_dir=$stacks_dir/populations_r80
 
mkdir -p ${out_dir}
 
log_file=$out_dir/populations.oe
 
populations -P $stacks_dir -O $out_dir -r 0.80 &> $log_file &
 
done

#!/bin/bash
 
# Extract the SNPs-per-locus distributions (they are reported in the log of populations).
 
# ----------
 
echo "Tallying the numbers..."
 
full_table=n_snps_per_locus.tsv
 
header=\'#par_set\tM\tn\tm\tn_snps\tn_loci\'
 
for M in 1 2 3 4 5 6 7 8 9 ;do
 
 n=$M
 
 m=3
 
 # Extract the numbers for this parameter combination.
 
 log_file=test/de-novo/stacks_M${M}/populations_r80/batch_1.populations.log
 
 sed -n \'/^#n_snps\tn_loci/,/^[^0-9]/ p\' $log_file | grep -E \'^[0-9]\' > $log_file.snps_per_loc
 
 # Cat the content of this file, prefixing each line with information on this
 
 # parameter combination.
 
 line_prefix="M$M-n$n-m$m\t$M\t$n\t$m\t"
 
 cat $log_file.snps_per_loc | sed -r "s/^/$line_prefix/"
 
done | sed "1i $header" > $full_table

SAMPLES, = glob_wildcards("01-clean-data/{sample}.fq.gz")
 
INDEX_DICT = {value: key for key, value in dict(enumerate(SAMPLES, start=1)).items()}
 
FQ_FILES = expand("01-clean-data/{sample}.fq.gz", sample=SAMPLES)
 
# de novo
 
MISMATCH = 4 # mismatch number of ustacks
 
DE_NOVO_LOCI = expand("03-stacks-analysis/de-novo/{sample}.snps.tsv.gz", sample=SAMPLES)
 
DE_NOVO_CATA = "03-stacks-analysis/de-novo/batch_1.catalog.snps.tsv.gz"
 
DE_NOVO_MATS = expand("03-stacks-analysis/de-nono/{sample}.matches.tsv.gz", sample=SAMPLES)
 
# ref-based
 
INDEX = "reference/index/bwa/gac"
 
BAM_FILES = expand("02-read-alignment/ref-based/{sample}.bam", sample=SAMPLES)
 
REF_BASED_LOCI = expand("03-stacks-analysis/ref-based/{sample}.snps.tsv.gz", sample=SAMPLES)
 
REF_BASED_CATA = "03-stacks-analysis/ref-based/batch_1.catalog.snps.tsv.gz"
 
REF_BASED_MATS = expand("03-stacks-analysis/ref-based/{sample}.matches.tsv.gz", sample=SAMPLES)
 
 
 
rule all:
 
 input: rules.de_novo.input, rules.ref_based.input
 
 
rule de_novo:
 
 input: DE_NOVO_LOCI, DE_NOVO_CATA,DE_NOVO_MATS
 
 
rule ref_based:
 
 input: REF_BASED_LOCI, REF_BASED_CATA, REF_BASED_MATS
 
 
# de novo data analysis
 
## The unique stacks program will take as input a set of short-read sequences
 
## and align them into exactly-matching stacks (or putative alleles).
 
rule ustacks:
 
 input: "01-clean-data/{sample}.fq.gz"
 
 threads: 8
 
 params:
 
 mismatch = MISMATCH,
 
 outdir = "03-stacks-analysis/de-novo",
 
 index = lambda wildcards: INDEX_DICT.get(wildcards.sample)
 
 output:
 
 "03-stacks-analysis/de-novo/{sample}.snps.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.tags.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.models.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.alleles.tsv.gz",
 
 log: "03-stacks-analysis/de-novo/{sample}.ustacks.oe"
 
 shell:"""
 
 mkdir -p {params.outdir}
 
 ustacks -p {threads} -M {params.mismatch} -m 3 \
 
 -f {input} -i {params.index} -o {params.outdir} &> {log}
 
 """
 
 
## choose sample for catalog building
 
rule choose_representative_samples:
 
 input: expand("03-stacks-analysis/de-novo/{sample}.ustacks.oe", sample=SAMPLES)
 
 output:
 
 "03-stacks-analysis/de-novo/per_sample_mean_coverage.tsv",
 
 "info/popmap.catalog.tsv"
 
 shell:"""
 
 for sample in {input};do
 
 name=${{sample##*/}}
 
 name=${{name%%.*}}
 
 sed -n "/Mean/p" ${{sample}} |\
 
 sed "s/.*Mean: \(.*\); Std.*/\\1/g" |\
 
 paste - - - |\
 
 sed "s/.*/${{name}}\\t&"
 
 done | sort -k2,2nr > {output[0]}
 
 head -n 50 {output[0]} | tail -n 40 | cut -f 1 | sed \'s/\([0-9a-zA-Z]*\)_.*/&\t\1/\' > {output[1]}
 
 """
 
 
rule de_novo_cstacks:
 
 input:
 
 "info/popmap.catalog.tsv",
 
 expand("03-stacks-analysis/de-novo/{sample}.snps.tsv.gz", sample=SAMPLES)
 
 params:
 
 stacks_path = "03-stacks-analysis/de-novo/",
 
 mismatch = MISMATCH
 
 threads: 10
 
 log:"03-stacks-analysis/de-novo/de_novo_cstacks.oe"
 
 output:
 
 "03-stacks-analysis/de-novo/batch_1.catalog.alleles.tsv.gz",
 
 "03-stacks-analysis/de-novo/batch_1.catalog.snps.tsv.gz",
 
 "03-stacks-analysis/de-novo/batch_1.catalog.tags.tsv.gz"
 
 shell:"""
 
 cstacks -p {threads} -P {params.stacks_path} -M {input[0]} \
 
 -n {params.mismatch} &> {log}
 
 """
 
 
rule de-novo_sstacks:
 
 input:
 
 "03-stacks-analysis/de-novo/{sample}.snps.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.tags.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.models.tsv.gz",
 
 "03-stacks-analysis/de-novo/{sample}.alleles.tsv.gz"
 
 params:
 
 catalog = "03-stacks-analysis/de-novo/batch_1",
 
 sample = lambda wildcards: "03-stacks-analysis/de-novo/"   wildcards.sample
 
 output:
 
 "03-stacks-analysis/de-novo/{sample}.matches.tsv.gz"
 
 log: "03-stacks-analysis/de-novo/{sample}.sstacks.oe"
 
 shell:"""
 
 sstacks -c {params.catalog} -s {params.sample} -o 03-stacks-analysis/de-novo &> {log}
 
 """
 
 
# reference-based data analysis
 
## read alignment with bwa-mem
 
rule bwa_mem:
 
 input: "01-clean-data/{sample}.fq.gz"
 
 params:
 
 index = INDEX,
 
 mismatch = "3",
 
 gap = "5,5"
 
 threads: 8
 
 output: "02-read-alignment/ref-based/{sample}.bam"
 
 shell:"""
 
 mkdir -p 02-read-alignment
 
 bwa mem -t {threads} -M -B {params.mismatch} -O {params.gap} {params.index} {input} \
 
 | samtools view -b > {output}
 
 """
 
 
## find variant loci
 
rule pstacks:
 
 input: "02-read-alignment/ref-based/{sample}.bam"
 
 params:
 
 outdir = "03-stacks-analysis/ref-based/",
 
 index = lambda wildcards: INDEX_DICT.get(wildcards.sample)
 
 threads: 8
 
 output:
 
 "03-stacks-analysis/ref-based/{sample}.snps.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.tags.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.models.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.alleles.tsv.gz"
 
 log: "03-stacks-analysis/ref-based/{sample}.pstacks.oe"
 
 shell:"""
 
 mkdir -p 03-stacks-analysis/ref-based
 
 pstacks -p {threads} -t bam -f {input} -i {params.index} -o {params.outdir} &> {log}
 
 """
 
 
## A catalog can be built from any set of samples processed by the ustacks or pstacks programs
 
rule ref_based_cstacks:
 
 input: "info/popmap.tsv",expand("03-stacks-analysis/ref-based/{sample}.snps.tsv.gz", sample=SAMPLES)
 
 threads: 10
 
 output:
 
 "03-stacks-analysis/ref-based/batch_1.catalog.alleles.tsv.gz",
 
 "03-stacks-analysis/ref-based/batch_1.catalog.snps.tsv.gz",
 
 "03-stacks-analysis/ref-based/batch_1.catalog.tags.tsv.gz"
 
 shell:
 
 "cstacks -p {threads} --aligned -P 03-stacks-analysis/ref-based/ -M {input[0]}"
 
 
## Sets of stacks, i.e. putative loci, constructed by the ustacks or pstacks programs
 
## can be searched against a catalog produced by cstacks.
 
rule ref_based_sstacks:
 
 input:
 
 "03-stacks-analysis/ref-based/{sample}.snps.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.tags.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.models.tsv.gz",
 
 "03-stacks-analysis/ref-based/{sample}.alleles.tsv.gz"
 
 params:
 
 catalog = "03-stacks-analysis/ref-based/batch_1",
 
 sample = lambda wildcards: "03-stacks-analysis/ref-based/"   wildcards.sample
 
 output:
 
 "03-stacks-analysis/ref-based/{sample}.matches.tsv.gz"
 
 log: "03-stacks-analysis/ref-based/{sample}.sstacks.oe"
 
 shell:"""
 
 sstacks --aligned -c {params.catalog} -s {params.sample} -o 03-stacks-analysis/ref-based &> {log}
 
 """

snakemake --cluster "qsub -V -cwd" -j 20 --local-cores 10 &

#!/usr/bin/env Rscript
 
 
snps_per_loc = read.delim(\'./n_snps_per_locus.tsv\')
 
# Keep only M==n, m==3
 
snps_per_loc = subset(snps_per_loc, M==n & m==3)
 
# Rename column 1
 
colnames(snps_per_loc)[1] = \'par_set\'
 
 
# Create a new data frame to contain the number of loci and polymorphic loci
 
d = snps_per_loc[,c(\'par_set\', \'M\', \'n\', \'m\')]
 
d = d[!duplicated(d),]
 
 
# Compute these numbers for each parameter set, using the par_set column as an ID
 
rownames(d) = d$par_set
 
for(p in rownames(d)) {
 
 s = subset(snps_per_loc, par_set == p)
 
 d[p,\'n_loci\'] = sum(s$n_loci)
 
 s2 = subset(s, n_snps > 0)
 
 d[p,\'n_loci_poly\'] = sum(s2$n_loci)
 
}
 
 
# Make sure the table is ordered
 
d = d[order(d$M),]
 
 
pdf(\'./n_loci_Mn.pdf\')
 
 
# Number of loci
 
# ==========
 
 
plot(NULL,
 
 xlim=range(d$M),
 
 ylim=range(c(0, d$n_loci)),
 
 xlab=\'M==n\',
 
 ylab=\'Number of loci\',
 
 main=\'Number of 80%-samples loci as M=n increases\',
 
 xaxt=\'n\',
 
 las=2
 
 )
 
abline(h=0:20*5000, lty=\'dotted\', col=\'grey50\')
 
axis(1, at=c(1,3,5,7,9))
 
legend(\'bottomright\', c(\'All loci\', \'Polymorphic loci\'), lty=c(\'solid\', \'dashed\'))
 
 
lines(d$M, d$n_loci)
 
points(d$M, d$n_loci, cex=0.5)
 
 
lines(d$M, d$n_loci_poly, lty=\'dashed\')
 
points(d$M, d$n_loci_poly, cex=0.5)
 
 
# Number of new loci at each step (slope of the previous)
 
# ==========
 
 
# Record the number of new loci at each parameter step
 
d$new_loci = d$n_loci - c(NA, d$n_loci)[1:nrow(d)]
 
d$new_loci_poly = d$n_loci_poly - c(NA, d$n_loci_poly)[1:nrow(d)]
 
 
# Record the step size
 
d$step_size = d$M - c(NA, d$M)[1:(nrow(d))]
 
 
plot(NULL,
 
 xlim=range(d$M),
 
 ylim=range(c(0, d$new_loci, d$new_loci_poly), na.rm=T),
 
 xlab=\'M==n\',
 
 ylab=\'Number of new loci / step_size (slope)\',
 
 main=\'Number of new 80%-samples loci as M=n increases\'
 
 )
 
abline(h=0, lty=\'dotted\', col=\'grey50\')
 
legend(\'topright\', c(\'All loci\', \'Polymorphic loci\'), lty=c(\'solid\', \'dashed\'))
 
 
lines(d$M, d$new_loci / d$step_size)
 
points(d$M, d$new_loci / d$step_size, cex=0.5)
 
 
lines(d$M, d$new_loci_poly / d$step_size, lty=\'dashed\')
 
points(d$M, d$new_loci_poly / d$step_size, cex=0.5)
 
 
null=dev.off()
 
plot_n_loci.R的代码
 
#!/usr/bin/env Rscript
 
 
d = read.delim(\'./n_snps_per_locus.tsv\')
 
# Keep only M==n, m==3.
 
d = subset(d, M==n & m==3)
 
# Make sure the table is ordered by number of snps.
 
d = d[order(d$n_snps),]
 
 
Mn_values = sort(unique(d$M))
 
 
# Write the counts in a matrix.
 
m = matrix(NA, nrow=length(Mn_values), ncol=max(d$n_snps) 1)
 
for(i in 1:nrow(d)) {
 
 m[d$M[i],d$n_snps[i] 1] = d$n_loci[i] # [n_snps 1] as column 1 is for loci with 0 SNPs
 
}
 
 
# Truncate the distributions.
 
max_n_snps = 10
 
m[,max_n_snps 2] = rowSums(m[,(max_n_snps 2):ncol(m)], na.rm=T)
 
m = m[,1:(max_n_snps 2)]
 
m = m / rowSums(m, na.rm=T)
 
 
# Draw the barplot.
 
pdf(\'n_snps_per_locus.pdf\')
 
 
clr = rev(heat.colors(length(Mn_values)))
 
 
barplot(m,
 
 beside=T, col=clr, las=1,
 
 names.arg=c(0:max_n_snps, paste(\'>\', max_n_snps, sep=\'\')),
 
 xlab=\'Number of SNPs\',
 
 ylab=\'Percentage of loci\',
 
 main=\'Distributions of the number of SNPs per locus\nfor a range of M==n values\'
 
 )
 
legend(\'topright\', legend=c(\'M==n\', Mn_values), fill=c(NA, clr))
 
 
null=dev.off()

# 项目根目录
min_samples=0.80min_maf=0.05max_obs_het=0.70populations -P 03-stacks-analysis/ref-based/ -r $min_samples --min_maf $min_maf \
--max_obs_het $max_obs_het --genepop &> populations.oe
# --genepop表示输出格式，可以选择--vcf等