遗传学基因组学重测序下游分析

受选择基因筛选plus—xpclr

2021-05-26  本文已影响0人  Morriyaty

之前讲过根据Fst和pi值筛选基因,加上一条根据xpclr筛选,最后根据三者的交集选出基因。下面讲一下xpclr的做法。

先下载xpclr
conda create -n xpclr 
conda activate xpclr
conda install xpclr
然后处理一下我们的VCF文件(项目示例中有8个群体)
1)提取只包含两个群体的VCF子文件
2)按染色体拆分
3)每个染色体VCF需要bgzip压缩并建立.tbi索引
4)将VCF和索引放入一个文件夹内
批量生成命令行
perl 0.xpclr.pl
在结果文件夹中:
1)log信息放入文件夹log中
2)结果文件放入文件夹res中
3)合并结果
perl 0.MergeRes.pl > xpclr.txt
提取前1%的染色体区域
cat xpclr.txt | grep -v "chrom" | sort -n -k 12 -r | sed 's/[0-9].[0-9]\+e-[0-9]\+/0/g' | sort -n -k 12 -r | head -n $1% | cut -f 1-3 | tr -s "\t" > top1.bed

后面的就很简单了,不提了。

######

##所有的 .log  放入log文件夹
##所有的 其他文件  放入res文件夹 

cp /data/01/user152/workspace/F/select-swip/gene.gff ./
cp /data/01/user152/workspace/F/select-swip/gene-pipei.list ./
cp /data/01/user152/workspace/F/select-swip/pipei.py ./

perl 0.MergeRes.pl > xpclr.txt

cat xpclr.txt | awk '{print$1,$2,$3,$12}' | sed 's/-[0-9.]/0/g' | awk  'IF$4>0{print$0}' | grep -v '0.0' | awk '{print$1"\t"$2"\t"$3"\t"$4}'| grep -v "chrom" | sort -t $'\t' -k 1,1 -k 2n,2 > xpclr_plot.txt

wc -l xpclr_plot.txt #5%

cat xpclr_plot.txt | sort -k 4 -nr | head -n 1740 | sort -t $'\t' -k 1,1 -k 2n,2> top5.txt

cat top5.txt | awk '{print$1"\t"$2"\t"$3}' | grep -v "chrom" > gene.bed

bedtools merge -i gene.bed > gene_merge.bed

bedtools  intersect -wo -F 0.1 -a gene_merge.bed -b gene.gff | awk -F "\t"  '{print $4"\t"$7"\t"$8"\t"$10"\t"$12}' | sort -u | awk '{print$5}' | awk -F ";" '{print$1}' | sed 's/ID=//g' > gene.list

bedtools  intersect -wo -F 0.1 -a gene_merge.bed -b gene.gff | awk '{print$1,$2,$3,$12}' > gene-xpclr.region

python3 pipei.py gene-pipei.list gene.list commongene.list
上一篇下一篇

猜你喜欢

热点阅读