07重测序流程2-遗传进化
一.准备数据(看4.染色体列表!!!)
1. 测序所在路径，一般有多个文件夹，每个样品一个文件夹
2. 基因组fasta文件，有时需要下载，有时客户提供，还有常见的或者我们做过的可以到 /share/nas6/pub/genome中找找
3. 基因组注释文件，gff3格式
$ perl /share/nas6/xul/program/chloroplast/mvista/genbank.to.gff.pl -i *.gbk -f *.fasta -o tmp.gff #gbk转gff
$ perl -lane '/Name=(.*;)/;print $_."Parent=$1ID=$1"' tmp.gff > genome.gff #修改名字
4. 染色体列表，txt文件，一列
存放序列名称，因为有的fasta文件中有很多序列，该文件用于定义某些图片中需要展示的序列。
一般是保留染色体级别的序列，可以保留所有染色体
如果该基因组都是scaffold，那么就保留几个比较长的序列名称，十个左右就行
$ ir.py -i *.fna -l | sort -k 1 -nr > len.txt
$ cat len.txt | head -n 15 | awk '{print $2}' > chrlist.txt && cat chrlist.txt # 根据情况取前多少行
5. 功能注释结果的整合文件
6. 分组group.txt 实际上用不到,因为会在配置文件那手动填写
组名 样品名
CQ CQ1
CS CS1
FJ FJ1
FJ FJ2
GX GX1
GX GX2
GZ GZ1
GZ GZ2
HB HB1
HN HN1
HN HN2
二.功能注释(有的可以跳过)
可放后台 如果已经做过了，可以跳过该步骤，nas6中的基因组一般都做过这个分析 /share/nas6/pub/genome
1. 获取转录本序列文件，fasta格式的核酸序列
gffread -x cds.fa -g genome.fa(参考基因组) in.gff3(gff3文件)
遇到错误就删删删
例子:
$ grep -n "rna-A564_p001" GCF_004118075.2_ASM411807v2_genomic.gff | wc -l
$ for i in {1..6};do echo $i;sed -i '649514d' GCF_004118075.2_ASM411807v2_genomic.gff;done
# cds要用id来命名，重测序这块分析可能不用管
2.获取最长转录本
有的基因可能存在多个转录本，所以需要过滤下，根据上述提取出来的基因名称判断哪些是同一个基因的，一般id后面有.1 .2这种，就属于一个基因的不同转录本。
具体还要看提取出来的格式。
程序：
# perl /share/nas6/zhouxy/functional_modules/cnvfmt_trans2longest/fasta2longest.pl -i cds.fa -o filter.cds.fa # 旧的程序不用
$ perl /share/nas6/zhouxy/functional_modules/cnvfmt_trans2longest/fasta2longest_with_gff3.pl -i cds.fa -g *.gff3 -o filter.cds.fa # 用这个
# perl /share/nas6/zhouxy/functional_modules/cnvfmt_trans2longest/fasta_trans_to_geneid_format_convert_ncbi.pl 格式转换 没用过
3. 填写功能注释的配置文件
# 植物 下面的nr库需要根据情况改成单子叶的或者双子叶的，默认双子叶
$ cp /share/nas6/pub/pipline/gene-function-annotation/v1.3/profile/plant_QAll.cfg ./
# 病毒
$ cp /share/nas6/pub/pipline/gene-function-annotation/v1.3/profile/VRL_QAll.cfg ./
# 真菌
$ cp /share/nas6/pub/pipline/gene-function-annotation/v1.3/profile/fungi_QAll.cfg ./
修改第一行，改成filter.cds.fa的绝对路径！
4. 运行主程序
$ nohup perl /share/nas6/pub/pipline/gene-function-annotation/v1.3/Gene_Func_Anno_Pipline.pl --nr --swissprot --cog --kog --kegg --pfam --GO --cfg *_QAll.cfg --od Basic_function &
程序运行完后需要用到的文件就是：Basic_function/Result/All_Database_annotation.xls 这个文件，用于填写下面的配置文件
三.主流程
/share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl
3.0 配置文件示例
3.0.1 有分组
/share/nas1/yuj/project/re-sequencing/GP-20220607-4434_GP-20220614-4481_1_20221124/ref_reseq_config.yaml2 # 有很多分组
多了Contrast和parameters
Contrast :
- CQ_vs_CS
- CQ_vs_FJ
- CS_vs_FJ
parameters:
windows : 100000
steps : 10000
Groups :
CQ : CQ1
CS : CS1
FJ : FJ1,FJ2
PSMC :
- CQ1
- CS1
- FJ1
- FJ2
3.0.2 无分组
/share/nas1/yuj/project/re-sequencing/GP-20230327-5799_20230519/e1345/ref_reseq_config.yaml # 无分组
Groups :
all : E1,E3,E4,E5
PSMC :
- E1
- E3
- E4
- E5
3.1 生成配置文件
perl /share/nas6/xul/program/reseq/create_profile_for_reseq_pop.pl # 有点问题
perl /share/nas1/yuj/program/reseq/create_profile_for_reseq_pop.pl
-i <in dir> contain all.gz file or fq or fastq
# 输入包含测序文件所在的路径（文件夹，不是具体文件路径），如 path/clean/*/
-s ref genome
# 输入参考基因组fasta文件
-g ref chr list
# 输入染色体列表
-c ref anno Gene_Anno/Result/All_Database_annotation.xls default
# 输入功能注释的结果Basic_function/Result/All_Database_annotation.xls，如果功能注释没做完，这个可以先不填。
-gff ref gff
# 输入gff文件
-Is insert_size,default 100-500
# 插入片段大小，一般不用改
-p Prefix default XM
# 物种名简拼，不改也行
# 手动添加分组 详见35个栀子分析 后续写程序改进一下
cd ../
perl /share/nas1/yuj/program/reseq/create_profile_for_reseq_pop.pl -i xxx/2.clean_data/*/ -s data/*.fna -g data/chrlist.txt -gff data/*.gff -c data/Basic_function/Result/All_Database_annotation.xls -p M8
# 直接生成的就是无分组
3.2 运行
主流程：
$ nohup perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis & # 公司
$ nohup perl /share/nas1/yuj/pipline/pop-pip/v1.1/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis & # 修改
gwas：/share/nas6/zhouxy/pipline/gwas_vcf_pip/current/gwas_ped_pip_v2.pl
3.3 报告
整理结果：
# perl /share/nas6/zhouxy/pipline/pop-pip/v1.1/script/resultdir/resultDir.pl -i analysis/ -o complete 原版有点问题,10.ssweep多打了字母和路径
$ perl /share/nas1/yuj/pipline/pop-pip/v1.1/script/resultdir/resultDir.pl -i analysis/ -o complete
报告配置：
$ cp /share/nas6/xul/project/ddrad/GP-20210729-3250_ningxiadaxue_211samples_niu_ddrad/pop_html.cfg . && realpath pop_html.cfg
生成报告:
# perl /share/nas6/zhouxy/pipline/pop-pip/v1.1/html_report/pop_report2html.pl -id complete/ -cfg pop_html.cfg 原版有点问题
$ perl /share/nas1/yuj/pipline/pop-pip/v1.1/html_report/pop_report2html.pl -id complete/ -cfg pop_html.cfg # 限制了表格的行数
五.FAQ
5.1 gatk_VQSR.failed.cmds
运行/share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.cmds_dir/gatk_VQSR.failed.cmds出现Exception in thread "main" java.lang.OutOfMemoryError: GC overhead limit exc
## 调大jvm的内存参数
--java-options "-Xmx64g -Xms64g"
5.2 gatk_Stat.failed.cmds
plot-vcfstats这一步 ---- gatk_Stat.failed.cmds中第一行命令的第二条命令
We are sorry, commands in file: [/share/nas1/yuj/project/re-sequencing/GP-20230517-6067_20230621/analysis/.cmds_dir/gatk_Stat.failed.cmds] failed. :-(
Parsing bcftools stats output: /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/variation_dir/variants_stat/samples.var.stat
Plotting graphs: python3 plot.py
Creating PDF: pdflatex summary.tex >plot-vcfstats.log 2>&1
The command exited with non-zero status, please consult the output of pdflatex: /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/variation_dir/variants_stat/samples.varplot-vcfstats.log
at /share/nas6/zhouxy/biosoft/bcftools/current/bin/plot-vcfstats line 111.
main::error("The command exited with non-zero status, please consult the o"...) called at /share/nas6/zhouxy/biosoft/bcftools/current/bin/plot-vcfstats line 2083
main::create_pdf(HASH(0xbb2d98)) called at /share/nas6/zhouxy/biosoft/bcftools/current/bin/plot-vcfstats line 72
Attaching package: ‘gplots’
The following object is masked from ‘package:stats’:
lowess
Warning message:
Removed 44825 row(s) containing missing values (geom_path).
Warning message:
Removed 44825 row(s) containing missing values (geom_path).
Error, cmd: ParaFly -c /share/nas1/yuj/project/re-sequencing/GP-20230517-6067_20230621/analysis/.cmds_dir/gatk_Stat.cmds -CPU 4 -shuffle -failed_cmds /share/nas1/yuj/project/re-sequencing/GP-20230517-6067_20230621/analysis/.cmds_dir/gatk_Stat.failed.cmds 2>tmp.105922.stderr died with ret 256 at /share/nas6/zhouxy/biosoft/perl/current/lib/site_perl/5.28.3/Pipeliner.pm line 102.
Pipeliner::run(Pipeliner=HASH(0x3619e68)) called at /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl line 1284
main::Pipline_parafly_commands("gatk_Stat.cmds", "gatk_Stat_analysis.ok", "gatk_Stat.failed.cmds", 4) called at /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl line 680
Perl exited with active threads:
1 running and unjoined
1 finished and unjoined
0 running and detached
## 具体的log
analysis/variation_dir/variants_stat/samples.var/plot-vcfstats.log
! LaTeX Error: File 'multirow.sty' not found.
5.2.1 pdflatex换成tectonic解决
1.画图
cd analysis/variation_dir/variants_stat/samples.var/ && python3 plot.py && /share/nas1/yuj/software/miniconda3/envs/bcftools/bin/tectonic summary.tex && cd -
2.运行gatk_Stat.failed.cmds中第一行命令的最后一部分
perl /share/nas1/yuj/pipline/pop-pip/v1.1/script/02variant/snp_ver.stat.pl -i analysis/variation_dir/variants_stat/samples.var.stat -o analysis/variation_dir/variants_stat/snp_var.stat.xls
3.制造一个gatk_Stat_analysis.ok接着跑
touch analysis/.flag_dir/gatk_Stat_analysis.ok && nohup perl /share/nas1/yuj/pipline/pop-pip/v1.1/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis &
安装tectonic:
conda create -n bcftools -y
conda install -c conda-forge Tectonic
conda activate bcftools
tectonic summary.tex # 生成summary.pdf
conda deactivate
pdflatex --version
tectonic --version
原因
plot-vcfstats脚本在bcftools安装目录misc文件夹下，这是一个perl脚本，会调用python3绘图模块Matplotlib。如果没有安装该模块，可以通过pip3命令pip3 install -U matplotlib进行安装。如果没有pip3，通过sudo apt-get install python3-pip命令安装。
plot-vcfstats生成pdf文件还需要pdf-latex，如果系统没有安装latex，通过sudo apt-get install texlive-full进行安装。
Quote
1."python - Dependencies plot-vcfstats in conda environment - Stack Overflow." https://stackoverflow.com/questions/69659849/dependencies-plot-vcfstats-in-conda-environment, 访问时间 1 Jan. 1970.
2."Plot-VCF issues · Issue #1697 · samtools/bcftools · GitHub." https://github.com/samtools/bcftools/issues/1697, 访问时间 1 Jan. 1970.
3."Impossible to include a PDF image with '.' in filename · Issue #187 · tectonic-typesetting/tectonic · GitHub." https://github.com/tectonic-typesetting/tectonic/issues/187, 访问时间 1 Jan. 1970.
5.3 vcf2plink.cmds
流程报错信息
structure: -- Too few valid variants for --indep-pairwise
!!! 与variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf有关
这一步第一行命令共3部分,产生最终的samples.plink.map
需要对samples.plink.map第一列编号进行处理
1. 生成id映射表 Create chromosome map
bcftools view -H analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf | cut -f 1 | uniq | awk '{count++; print $0"\t"count}' > analysis/PTS_analysis/admixture_dir/chrom-map.txt
chrom-map.txt样式:
LUGG01000001.1 1
LUGG01000002.1 2
LUGG01000003.1 3
LUGG01000004.1 4
LUGG01000005.1 5
2. Make ped file using this chromosome map
vcftools --vcf analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf --plink --chrom-map analysis/PTS_analysis/admixture_dir/chrom-map.txt --out analysis/PTS_analysis/admixture_dir/samples.plink
# 运行后生成 samples.plink.map
3.运行vcf2plink.cmds中的第二行命令
## 如果染色体条数小于23，直接运行第二行，不需要更改
sed -n '2p' analysis/.cmds_dir/vcf2plink.cmds | sh
## 如果染色体条数超过23，第二行每一步都要加上 --allow-extra-chr
plink --allow-extra-chr
4.
touch analysis/.flag_dir/vcf2plink.ok
5.继续运行主流程
perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
！！！这下面的不用了哈！！！
0.直接运行第一条命令的第一部分,生成初始的samples.plink.map
$ vcftools --vcf analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf --plink --out analysis/PTS_analysis/admixture_dir/samples.plink
1.下面的命令用来显示map文件里染色体有哪些
如果没有按自然顺序排序,结果会比较奇怪
$ cd analysis/PTS_analysis/admixture_dir
$ time (a="start" && for i in `awk '{print $2}' samples.plink.map | awk -F ':' '{print $1}'`;do if [ $a != $i ];then echo $i ;a=$i;fi;done | sort -k 1,1n > all.id) #awk可以
$ time (a="start" && for i in `cut -f 2 -d
## 5.4 vcffiles2plink.cmds 原因同5.3
```shell
0.生成id映射表
bcftools view -H analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf | cut -f 1 | uniq | awk '{count++; print $0"\t"count}' > analysis/PTS_analysis/pca_dir/chrom-map.txt
chrom-map.txt格式如下:第一列id 第二列替换后的id(用数字表示)
LUGG01000001.1 1
LUGG01000002.1 2
LUGG01000003.1 3
LUGG01000004.1 4
LUGG01000005.1 5
LUGG01000006.1 6
LUGG01000007.1 7
LUGG01000009.1 9
LUGG01000010.1 10
...
LUGG01000127.1 46
1.直接生成samples.plink.map
vcftools --vcf analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf --plink --chrom-map analysis/PTS_analysis/pca_dir/chrom-map.txt --out analysis/PTS_analysis/pca_dir/samples.plink
2.运行vcffiles2plink.cmds的后三行命令
小于23,直接运行后三行,不需要更改
sed -n '4p' /share/nas1/yuj/project/re-sequencing/GP-20230912-6827_20230914/analysis/.cmds_dir/vcffiles2plink.cmds | sh && sed -n '5p' /share/nas1/yuj/project/re-sequencing/GP-20230912-6827_20230914/analysis/.cmds_dir/vcffiles2plink.cmds | sh && sed -n '6p' /share/nas1/yuj/project/re-sequencing/GP-20230912-6827_20230914/analysis/.cmds_dir/vcffiles2plink.cmds | sh
如果染色体条数超过23,plink每一步都加上 --autosome-num 46 # 这里46指该物种的染色体条数
3.
touch analysis/.flag_dir/plink2mkbed.ok
4.继续运行主流程
perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
5.当染色体条数多于23时,接下来的报错
Pipline_sh_commands("gcta_pca.cmds","gcta_pca.ok");
sh analysis/.cmds_dir/gcta_pca.cmds
涉及到gcta的步骤都要加上 --autosome-num 46
touch analysis/.flag_dir/gcta_pca.ok
5.5 fstpi无结果(直接跳过即可)
analysis/.cmds_dir/selective_sweep_filter.failed.cmds出错
分别只有一个样品的三组组间互相对比无结果
touch /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.flag_dir/selective_sweep_filter.ok
perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
selective_sweep_region.failed.cmds 同上
$ touch /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.flag_dir/selective_sweep_region.ok
$ perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
selective_sweep_info.failed.cmds 同上
$ touch /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.flag_dir/selective_sweep_info.ok && perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
selective_sweep_anno.failed.cmds 同上
$ touch /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.flag_dir/selective_sweep_anno.ok && perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
5.6 step 14 vcf2plink2_for_gwas_input.cmds 同5.3
# vcf2plink2_for_gwas_input.cmds出错
# 命令1出错,同5.3,但输入文件变成以下文件
# input: analysis/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf
# out: analysis/GWAS_input/samples.plink.map
# 命令2也会出错,要改map中第一列
0.删除vcf中非染色体的序列
mv analysis/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf analysis/variation_dir/variants_anno_dir/backup_samples.pop.snp.anno.result.vcf # 备份
grep "##contig=<ID=" analysis/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf | awk -F "##contig=<ID=" '{print $2}' | awk -F "," '{print $1}' > id_contig.txt # 修改,留下需要删除的id
a='' && for i in `cat id_contig.txt`;do b=`echo "$a""\|""$i"`;a=$b;done && echo $b # 生成sed '/所需参数/d' file命令中的参数,前面会多一个 \| ,记得删除
sed '/参数/d' analysis/variation_dir/variants_anno_dir/backup_samples.pop.snp.anno.result.vcf > analysis/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf
# 例子 sed '/utg11958_G9\|utg18436_G9\|utg19729_G9\|utg7422_G9\|utg8976_G9\|utg9702_G9\|utg9928_G9\|utg1208_G1\|utg18778_G1\|utg3548_G1\|utg6849_G1\|utg7335_G1\|utg7888_G1\|utg9700_G1\|utg10062_G4\|utg13405_G4\|utg14560_G4\|utg16235_G4\|utg18525_G4\|utg18682_G4\|utg18967_G4\|utg2945_G4\|utg33336_G4\|utg42970_G4\|utg6726_G4\|utg80_G4\|utg1669_G2\|utg18871_G2\|utg6786_G2\|utg9385_G2\|utg9405_G2\|utg14488_G11\|utg34055_G11\|utg10629_G7\|utg10757_G7\|utg43111_G7\|utg6563_G7\|utg7405_G7\|utg7585_G7\|utg9635_G7\|utg1693_G8\|utg26320_G8\|utg28933_G8\|utg9041_G8\|utg10452_G5\|utg12633_G5\|utg13932_G5\|utg16492_G5\|utg23470_G5\|utg38007_G5\|utg4968_G5\|utg7033_G5\|utg7599_G5\|utg8857_G5\|utg1478_G3\|utg3749_G3\|utg6883_G3\|utg8759_G3\|utg11280_G6\|utg13378_G6\|utg20794_G6\|utg11850_G10\|utg14608_G10\|utg3856_G10\|utg8859_G10\|utg8952_G10\|utg9527_G10\|utg9722_G10\|utg12567\|utg18663/d' analysis/variation_dir/variants_anno_dir/backup_samples.pop.snp.anno.result.vcf > analysis/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf
1.重新运行第14步
sh /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.cmds_dir/vcf2plink2_for_gwas_input.cmds
# 也就是运行主流程 perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
接下来会有新的出错
2.命令2出错,对samples.plink.map进行手动替换
# map文件第一列以染色体编号做序号
# 使用for循环批量处理
$ cd /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/GWAS_input
$ for i in {1..11};do echo sed -i "'s/0\tGardenia$i:/$i\tGardenia$i:/g'" samples.plink.map;done > sed.sh # 显示ok
$ sh sed.sh
3.对samples.plink.map进行自然排序 -n
$ mv samples.plink.map tmp.samples.plink.map
$ cat tmp.samples.plink.map | sort -k 1 -n > samples.plink.map
4.运行第二条命令
5.
$ touch /share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/.flag_dir/vcf2plink2_for_gwas_input.ok
6.继续运行主流程
$ perl /share/nas6/zhouxy/pipline/pop-pip/current/pop_pip_v5.pl -i ref_reseq_config.yaml -o analysis
5.7 PSMC有效群体大小
另外一个程序
perl /share/nas1/yuj/pipline/genetic_diversity_pip/v1.0/script/06psmc/vcf2psmcGroup.pl -i analysis/vcf_filter/samples.pop.snp.recode.vcf -o analysis/psmc_analysis -g analysis/config_dir/group.list &
六.GWAS
perl /share/nas6/zhouxy/pipline/gwas_vcf_pip/current/gwas_ped_pip_v2.pl -i ddrad_config.yaml -o analysis
\t' samples.plink.map | cut -f 1 -d ':' `;do if [ $a != $i ];then echo
2.对samples.plink.map进行手动替换,map文件第一列以染色体编号做序号 不要用这种方法了,多次遍历浪费时间
#以/share/nas1/yuj/project/GP-20220607-4434_GP-20220614-4481_1_20221124/analysis/PTS_analysis/admixture_dir/samples.plink.map为例,以下三种命令都能替换成功:
perl -pi -e "s/0\tGardenia10/10\tGardenia10/g" samples.plink.map # ok
perl -pi -e 's/0\tGardenia9/9\tGardenia9/g' samples.plink.map # ok
sed -i 's/0\tGardenia8/8\tGardenia8/g' samples.plink.map # ok
#使用其中一种方式,for循环批量处理,例如:
for i in {1..11};do echo sed -i "'s/0\tGardenia
$ for i in {148..201};do echo sed -i.bak "'s/0\tNW_006711这是替换前的内容 #
$ sh sed.sh
3.对samples.plink.map进行自然排序 -n
$ mv samples.plink.map tmp.samples.plink.map
$ cat tmp.samples.plink.map | sort -k 1 -n > samples.plink.map
## 5.4 vcffiles2plink.cmds 原因同5.3
{{CODE_BLOCK_16}}
[vcf转化为plink报错](../生信软件/vcf转化为plink报错.md)
## 5.5 fstpi无结果(直接跳过即可)
{{CODE_BLOCK_17}}
selective_sweep_region.failed.cmds 同上
{{CODE_BLOCK_18}}
selective_sweep_info.failed.cmds 同上
{{CODE_BLOCK_19}}
selective_sweep_anno.failed.cmds 同上
{{CODE_BLOCK_20}}
## 5.6 step 14 vcf2plink2_for_gwas_input.cmds 同5.3
{{CODE_BLOCK_21}}
## 5.7 PSMC有效群体大小
{{CODE_BLOCK_22}}
# 六.GWAS
{{CODE_BLOCK_23}}