生信作业20220110

# 1. 从拟南芥的注释结果中提取GASA基因家族的seqid 
cat Athaliana_167_TAIR10.annotation_info.txt |grep GASA| awk -F "\t" '{print $3}' > At.GASA.list
# 2. 根据seqid，使用seqtk工具，提取相应的蛋白fasta序列
seqtk subseq Athaliana_167_TAIR10.protein.fa At.GASA.list > At.GASA.fa
# 3. 比对拟南芥的GASA蛋白序列
muscle -in At.GASA.fa  -out At.GASA.align.fa
# 4. 根据拟南芥的GASA基因家族比对结果，构建hmm模型
hmmbuild At.GASA.hmm At.GASA.align.fa
# 5. 在大豆的所有蛋白序列中，进行hmm检索
hmmsearch --tblout Gmax.GASA.list At.GASA.hmm Gmax_275_Wm82.a2.v1.protein.fa
# 6. 获取大豆的第一次GASA基因家族检索结果
awk '/^G/{print $1}' Gmax.GASA.list > Gmax.GASA.1st.list
seqtk subseq Gmax_275_Wm82.a2.v1.protein.fa Gmax.GASA.1st.list > Gmax.GASA.1st.fa
# 7. 根据第一次检索结果进行进一步建模，扩大检索范围，获取大豆的第二次GASA基因家族检索结果
muscle -in Gmax.GASA.1st.fa -out Gmax.GASA.1st.align.fa
hmmbuild Gmax.GASA.1st.hmm Gmax.GASA.1st.align.fa
hmmsearch --tblout Gmax.GASA.list.2 Gmax.GASA.1st.hmm Gmax_275_Wm82.a2.v1.protein.fa
awk '/^G/{print $1}' Gmax.GASA.list.2 > Gmax.GASA.2nd.list
seqtk subseq Gmax_275_Wm82.a2.v1.protein.fa Gmax.GASA.2nd.list > Gmax.GASA.2nd.fa
# 8. 合并拟南芥和大豆的GASA蛋白序列
cat Gmax.GASA.2nd.fa At.GASA.fa > At.Gmax.GASA.fa
muscle -in At.Gmax.GASA.fa -OUT At.Gmax.GASA.align.fa
trimal -in At.Gmax.GASA.align.fa -out At.Gmax.GASA.align.trimmed.fa -automated1

# 由于建树时，显示距离太大，无法建树，因而对上游结果进行检查
mkdir 2nd
cp *2nd
cd 2nd
awk '!/^>/{print length($0)}' Gmax.GASA.1st.fa|sort |uniq
# 发现有条序列过于长，应该去除该序列
awk 'NF%2==0{ID=$0}NF%2==1{len=length($0);if(len>250){print ID}}' Gmax.GASA.1st.fa
echo "Glyma.09G123000.1.p" > tmp.lst
# 删除该序列
sed -i '/Glyma.09G123000.1.p/d' Gmax.GASA.2nd.list
seqtk subseq Gmax_275_Wm82.a2.v1.protein.fa Gmax.GASA.1st.list > Gmax.GASA.1st.fa
muscle -in Gmax.GASA.1st.fa -out Gmax.GASA.1st.align.fa
hmmbuild Gmax.GASA.1st.hmm Gmax.GASA.1st.align.fa
hmmsearch --tblout Gmax.GASA.list.2 Gmax.GASA.1st.hmm Gmax_275_Wm82.a2.v1.protein.fa
awk '/^G/{print $1}' Gmax.GASA.list.2 > Gmax.GASA.2nd.list
seqtk subseq Gmax_275_Wm82.a2.v1.protein.fa Gmax.GASA.2nd.list > Gmax.GASA.2nd.fa
awk 'NF%2==0{ID=$0}NF%2==1{len=length($0);if(len>250){print ID}}' Gmax.GASA.2nd.fa
echo "Glyma.09G123000.1.p" >> tmp.lst
sed -i '/Glyma.02G145700.1.p/d' Gmax.GASA.2nd.list
seqtk subseq Gmax_275_Wm82.a2.v1.protein.fa Gmax.GASA.2nd.list > Gmax.GASA.2nd.fa
cat Gmax.GASA.2nd.fa At.GASA.fa > At.Gmax.GASA.fa
awk 'NF%2==0{ID=$0}NF%2==1{len=length($0);if(len>250){print ID}}' At.Gmax.GASA.fa
muscle -in At.Gmax.GASA.fa -OUT At.Gmax.GASA.align.fa
trimal -in At.Gmax.GASA.align.fa -out At.Gmax.GASA.align.trimmed.fa -automated1

awk '$1 ~ /^G/{print $1"\tT1";} $1 ~ /^A/{print $1"\tT2"}' evol.txt
awk '/AT2G14900.1/ { getline; print length($0) }' At.GASA.fa