仅作个人笔记使用意思就是只考虑我能不能看懂
参考
https://www.jianshu.com/p/3822d736234c
https://www.jianshu.com/p/ae320da46b29
先上这个数据库下载模式物种的目标基因或者基因家族
https://www.ebi.ac.uk/interpro/
1
2
3
当然还能用这个下载下来的pep文件去blast比对，这里先介绍hmmer
mafft --auto --clustalout OR.review.pep > OR.review.pep.clustal
上这个网站把mafft的结果转为stockholm格式
http://sequenceconversion.bugaco.com/converter/biology/sequences/fasta_to_phylip.php
4
hmmbuild OR.review.hmm input.stockholm
hmmsearch OR.review.hmm in.pep > hmmer.out
cat hmmer.out |grep -v "#" | awk '{if($4 < 1e-5 && $5 > 90) print $9}' | sort | uniq | grep -v "+" > hmmer.filter.out
blast比对
makeblastdb -in OR.review.pep -dbtype prot -out OR.review
blastp -num_threads 5 -db OR.review -query input.pep -outfmt 7 -seg yes > tmp
grep -v '#' tmp > .blastp.out
#筛选identity大于75%的序列
cat blastp.out |awk '$3>75' |cut -f1 |sort -u > blastp_result_id.list
将上述两种方法得到gene id合并取交集，找出两种方法共有的基因家族成员，使结果更可信
comm -12 blastp_result_id.list hmmer.filter.ou t> common.list
后面需要计算这些基因上的Fst
用的是pixy计算, 画1k的窗口
sort -k1,1 -k2,2n pixy.out > tmp
/data/01/user164/software/bedGraphToBigWig tmp chr.length car.vs.jud.Chr.Fst.bw
/data/01/user164/software/bigWigAverageOverBed car.vs.jud.Chr.Fst.bw OR.gene.bed OR.gene.Fst
awk '{if($3>0)print}' OR.gene.Fst | awk '{print "OR""\t"$6}' > OR.plot
用二项分布检验来做统计显著性
#全基因组随机抽3000个窗口或者更多窗口，计算平均值
#重复上述步骤1000次
#用GPT写了一个perl脚本
#GPT写的perl脚本比python脚本好使是怎么回事？？
perl random.pl tmp 3000 1000 | awk -F ": " '{print $2}' | sed 's/Mean = //g' | awk '{print "random""\t"$1}' > random.plot
vim random.pl
输入文件只有一列, 也就是每个窗口的Fst值
#!/usr/bin/perl
use strict;
use warnings;
# 从命令行参数中获取文件名、抽样行数和迭代次数
my $file_name = $ARGV[0];
my $num_samples = $ARGV[1];
my $num_iterations = $ARGV[2];
# 读取文件并将每行数据存储在数组中
open(my $file, '<', $file_name) or die "Cannot open file: $!";
my @lines = <$file>;
close($file);
# 删除每行末尾的换行符，并将字符串转换为数字
chomp(@lines);
my @data = map { $_ + 0 } @lines;
# 获取数据长度
my $data_len = scalar @data;
# 迭代执行第一步和第二步
for (my $i = 0; $i < $num_iterations; $i++) {
# 随机抽取指定行数的数据索引
my @sample_indices = pick_random_indices($num_samples, $data_len);
# 通过索引获取抽样数据
my @sample = @data[@sample_indices];
# 计算平均值
my $mean = calculate_mean(@sample);
# 打印结果
print "Iteration " . ($i + 1) . ": Mean = $mean\n";
}
# 选择指定数量的随机索引
sub pick_random_indices {
my ($count, $max) = @_;
my %picked;
my @indices;
while (scalar @indices < $count) {
my $index = int(rand($max));
push @indices, $index unless $picked{$index};
$picked{$index} = 1;
}
return @indices;
}
# 计算平均值
sub calculate_mean {
my @values = @_;
my $sum = 0;
foreach my $value (@values) {
$sum += $value;
}
return $sum / scalar(@values);
}
