基因家族分析 | 番茄Nramp基因家族分析（一）

以下分析完全参照lxmic的简书内容（基于全基因组的基因家族分析），侵删~

1.准备数据

①下载番茄全基因组CDS序列、protein序列、gff注释文件、基因组fa文件

ITAG3.2_CDS.fa
ITAG3.2_proteins.fa
ITAG3.2_gene_models.gff
S_lycopersicum_chromosomes.3.00.fa

②选定目标基因家族
目标基因家族： The natural resistance-associated macrophage protein (NRAMP)

③下载NRAMP的hmm模型 （https://pfam.xfam.org/）

Nramp.hmm   #PF01566

2.利用hmmsearch初步筛选基因家族成员

hmmsearch Nramp.hmm  Solanum_lycopersicum.SL2.50.pep.all.fa> out.txt
less -S  out.txt  
#下图结果可以看到在番茄蛋白序列中找到了9个属于Nramp基因家族的成员

番茄Nramp基因家族成员

3.批量获取家族成员信息

# 截取id号
vim out.txt
# 获取id号所在的行号，然后再用sed命令截取行，再用grep命令将id号匹配并重定向。
在vim命令模式下，输入“：set nu”
# sed命令截取，并用管道符直接输入给grep，匹配重定向到id文件
sed -n '17,26p' out.txt | grep -o "Sol.*\.1" > id.txt
# 利用samtools工具来进行序列提取# 首先建立索引文件
/data/zhaoqingyuan/biosoft/samtools-1.3.1/samtools faidx Solanum_lycopersicum.SL2.50.pep.all.fa    (最近镜像出现问题不能直接调用samtools)
# 再将id好作为输入，之后在重定向
xargs /data/zhaoqingyuan/biosoft/samtools-1.3.1/samtools faidx Solanum_lycopersicum.SL2.50.pep.all.fa < id.txt > nramp_protein_1.fasta &
less nramp_protein_1.fasta
# 得到的序列文件是含有回车符的，我利用一个perl单行命令将fasta格式的多行序列变成单行的fasta格式序列，链接：http://www.biotrainee.com/thread-291-1-1.html
erl -pe '/^>/ ? print "\n" : chomp' nramp_protein_1.fasta | tail -n +2 > nramp_protein.fasta
# 最后在samrt网站确认是否是该家族成员，进行最后的鉴定。

4.CDS序列获取

# Nramp_num文件含有id号
/data/zhaoqingyuan/biosoft/samtools-1.3.1/samtools faidx Solanum_lycopersicum.SL2.50.cds.all.fa
xargs /data/zhaoqingyuan/biosoft/samtools-1.3.1/samtools faidx Solanum_lycopersicum.SL2.50.cds.all.fa < id.txt > Nramp_cds_1.fasta
less Nramp_cds_1.fasta
# 发现CDS序列有回车，分为很多行，用Perl单行命令将其变成一行
perl -pe '/^>/ ? print "\n" : chomp' Nramp_cds_1.fasta | tail -n +2 > Nramp_cds.fasta
less Nramp_cds.fasta

5.Genomic DNA序列获取

# 根据geneid批量获取基因的染色体位置信息文件
grep -f id.txt ITAG3.2_gene_models.gff | cut -f 1,4,5,9 | grep "ID=mRNA" | sed 's/;Name.*//g' | sed 's/ID.*://g' > gene.bed
# 利用bedtools工具，将基因序列提取出来并重定向到新文件
/data/zhaoqingyuan/biosoft/samtools-1.3.1/samtools faidx Solanum_lycopersicum.SL2.50.dna.toplevel.fa
bedtools getfasta -fi S_lycopersicum_chromosomes.3.00.fa -bed gene.bed -name > gene_seq.fasta    (bedtools直接安装在环境中)

---

---

通过以上分析我们得到了以下数据：

id.txt
Nramp_cds.fasta
nramp_protein.fasta
gene.bed
gene_seq.fasta

接下来利用获得的这些数据进行可视化。

---

参考：
01 基于全基因组的基因家族分析（1）：数据准备
02 基于全基因组的基因家族分析（2）：SlNRAMP家族基因成员鉴定
03 基于全基因组的基因家族分析（3）：SlNRAMP家族基因CDS和Genomic DNA序列获取