仅作个人笔记 - 单基因家族分析

仅作个人笔记使用意思就是只考虑我能不能看懂

参考
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);
}