用find提取所有Dre和Ome的文件
find ./ -type f -name “gene.*” -exec cat {} + >gene.txt
find ./ -type f -name “pro.*” -exec cat {} + >pro.txt
剔除无用信息，杂乱id
#删除下载到的文件中注释行其他无用的信息
awk '/^>/ {split($2, arr, "["); print ">" arr[1]; next} {print}' Ome_pro.txt > Ome_proID.txt
# /^>/：匹配以 > 开头的行。
# split($2, arr, "[")：将第二列的内容按 [ 分割，取第一个部分，这样我们就提取了 LOC ID。
# print ">" arr[1]：只打印 LOC ID，并在前面加上 >。
# next：跳过其他的处理，只对以 > 开头的行执行上述操作。
# {print}：输出其他蛋白行。合并
把下载的比对文件合并为一个fatsa用于比对
cat ../../Dre_CYP家族成员/Dre_pro.txt ../../Ome_CYp家族成员/Ome_pro.txt >Dre_Ome_pro.fast
把湛江的（转录本后缀）删掉，处理一下id
sed 's/(.*)//g' needblast.fasta > cleaned_blast.fasta
建库
hmmsearch --cut_tc --domtblout ./zj/zj.PF.domtblout -o ./zj/zj.hmmout PF00067_seed.hmm ~/genome/zj_rmTE.aa
makeblastdb -in ~/genome/zj_rmTE.aa -dbtype prot -out zj_p
hmm模型构建，提取E小于0.05的行
hmmsearch --cut_tc --domtblout ./zj/zj.PF.domtblout -o ./zj/zj.hmmout PF00067.hmm ~/genome/zj_rmTE.aa
hmmsearch --domtblout ./zj/zj.PF.domtblout -o ./zj/zj.hmmout PF00067.hmm ~/genome/zj_rmTE.aa
awk '$7< 0.05 && $1 !~ /^#/ {print $0}' ./zj/zj.PF.domtblout > ./zj/zj.domtblout.filter
blast比对
blastp -query needblast.fasta -db zj_p -evalue 1e-5 -outfmt '6 std qlen slen' -out zj_p_id.blastout &
这里的needblast将弓背P450和Dre，Ome一起合并。
过滤结果
awk ' $3 > 30 ' zj_p_id.blastout > zj.blastout.filter.id
取并集，去重
cat ../../CYP-PF模型/zj/zj.filter.id ../1.blast/zj.blastout.filter.id | sort -u > zj_p_jiaoji.geneID
提取最长转录本
perl ~/Downloads/CYP/script/gff_filter_longest.pl zj_rmTE.gff3 zj_longest.txt zj_gff3
#用脚本提，第一个gff3是输入文件，后面txt是最长转录本id，gff3是过滤后gff3,此时你获得的txt是整个基因组的最长转录本，需要你使用你得到的基因家族成员的基因名去筛选提取出对应的最长转录本。
#找对应mapping,筛选出是你需要的转录本。
#R语言版本######################
#这里用R，这里的mapping.txt不是下文的，而是单独的一个只有geneid的txt，用来提取最长转录本的
#library(dplyr)
#result <- longest_need %>%
# filter(V1 %in% mapping$V1) %>%
# select(V2)
#
#write.table(result,"~/Downloads/CYP/sy/3.sort/longestid.txt",row.names = FALSE,,col.names = FALSE,quote = FALSE)
#################################
#linux语言版本###################
#也可以使用linux提最长转录本
#直接将你的基因家族成员ID复制为一列，创建一个mapping.txt
#grep -Ff mapping.txt zj_longest.txt >longestid.txt
#################################
这样就从zj_longest.txt中得到了呢基因对应的最长转录本,两个得到的结果是一致的，如果在服务器或MacOS上使用linux更方便，在个人windows计算机上使用R会更方便。
#使用seqkit提取最长转录本对应的蛋白序列
seqkit grep -f zj_longestid.txt zj_rmTM.aa >zj_pro.fasta
#得到最长转录本cds的蛋白序列提取到excel做表
转录本和基因名提取出来，做一个mapping.txt，这个mapping.txt第一列为基因ID，第二列为蛋白成员名，并做一个基因家族成员的统计表。
将geneID改名为蛋白名
#有些蛋白文件是先有.t1后缀的注释
sed -E 's/(>g[0-9]+)\.t[0-9]+/\1/' sy_pro.fasta > cleaned_sy.fasta
#这样出来的就没注释了，可以用mapping.txt后续替换成蛋白名称
awk 'NR==FNR {map[$1] = $2; next} /^>/ {sub($1, ">" map[substr($1, 2)]);}1' mapping.txt zj.fa > zj_gene_id.fa
# NR==FNR {map[$1] = $2; next}：当处理 mapping.txt 时，建立一个映射表，将第一列的基因ID映射到第二列的蛋白名称。
# /^>/ {sub($1, ">" map[substr($1, 2)]);}：对于 protein_sequences.fa 中的每一行，如果行首是 >，则根据映射表将基因ID替换为对应的蛋白名称。
# 1：表示打印替换后的行。awk '/^>/ {if (NR!=1) printf("\n"); printf("%s\n", $0); next;} {printf("%s", $0);} END {printf("\n");}' zj_P450_gene_副本.txt > output_file.txt
#去除蛋白序列文件中的换行符awk '/^>/ {if (seq) {print id "\t" seq;}id = substr($1, 2); seq = ""; }/^[^>]/ {seq = seq $0; }END {print id "\t" seq; }'"P450蛋白序列，去换行符.txt" > output_file.txt
#把蛋白序列从行变为列将蛋白序列放入表里，对应关系，计算长度。
-等电点分析
Rscript ~/Documents/CYP/script/fasta_pI_mw.R input output 合并需要作图的样本
多物种的时候需要合并为一个文件fasta文件，但物种数量过多的时候，mega分析会闪退崩溃。
cat P450.fa zj_gene_id.fa > mega.fasta
muscle处理，对齐序列
muscle -in mega.fa -out muscle_mega.fa
用mega出图
使用MODELS选项，进入find best DNA/PROTEIN
Models(ML)，使用ML法寻找最佳模型，并构建模型，重复次数1000次
当序列过大，数量过多时mega会崩溃，使用服务器软件进行拟合出图。
序列过大如果不用MEGA作图，使用服务器作图的话
多序列处理muscle -in all_Ome.fasta -out all_Ome.pep.mfa
NJ建树treebest nj all_Ome.pep.mfa -W -t jtt -b 100 > all_Ome.treebest.out
建树处理sed 's/:\([0-9.]\+\)\[&&NHX:B=\([0-9]\+\)\]/\2:\1/' \ all_Ome.treebest.out | awk '{printf $0}' > all_Ome.treebest.nwk
iqtrre自动计算模型建树
iqtree -s muscle.mfa -bb 1000 -m MFP -nt 20 -pre all_muscle
基因结构预测
..script/gff_filter_bymRNAID.pl genome/dmy_deletion_zj_chr.gff3 /CYP/nodmy/3.sort/longestid.txt tmp nodmy.gff3
提取转录本
去除无用碱基
's/\t\S\+UTR\t/\tUTR\t/i' nodmy.gff3 > nodmy.UTR.gff3
点击该网站进行预测
基因簇绘图
# 严格匹配一整列基因ID，要导入目标物种的gff3文件
使用上面已经提取过最长转录本的gff3文件
sed -E 's/ID=g([0-9]+)\.t[0-9]+.*/g\1/g' zj_id.gff3 |awk '$3 == "mRNA"'
#过去掉后缀，只选用mRNA列 将转录本中有对应的基因id提取出来，用gggenes画图
# 1. 将 V9 列重命名为 V2
matches$V2 <- matches$V9
# 2. 删除 V9 列中的 "ID=" 和 ";"
matches$V2 <- gsub("ID=", "", matches$V2) # 删除 "ID="
matches$V2 <- gsub(";", "", matches$V2) # 删除 ";"
# 3. 删除原有的 V2, V3, V6, V8 列
matches <- matches[, !colnames(matches) %in% c( "V3", "V6", "V8","V9")]
# 将 V7 列中的 "-" 替换为 "reverse"，"+" 替换为 "forward"
matches$V7 <- gsub("-", "reverse", matches$V7) # 替换 "-"
matches$V7 <- gsub("\\+", "forward", matches$V7) # 替换 "+"
library(gggenes)
ggplot(example_genes, aes(xmin = start, xmax = end, y = molecule, fill = gene)) + #x轴，y轴
geom_gene_arrow() +
facet_wrap(~ molecule, scales = "free", ncol = 1) +
scale_fill_brewer(palette = "Set3") +
theme_genes()转为gtf
gffread -T -o zj_id.gtf zj_id.gff3
提取上游2000bp
seqkit subseq --gtf dmy_id.gtf --feature transcript --up-stream 2000 --id-ncbi --only-flank --gtf-tag transcript_id ~/Downloads/genome/filter.new.genome.fa > dmy.script.fasta点击此处进入顺式作用元件预测
#首先有过滤后gff3，以及蛋白文件。这个方法处理后，得到的是只有目标基因家族成员的共线性分析
#建库
`makeblastdb -in ../nodmy_pro.fasta -dbtype prot`
`blastp -query ../nodmy_pro.fasta -db nodmy_pro.fasta -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out nodmy.blast`
nohup blastp -query Ocu_all.aa -db Ocu -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out Ocu.blast -num_threads 30 &
#gff3信息过滤,变为gff
`awk '$3 =="mRNA" ' ../zj_filter.gff3 | awk -F ";" '{print $1}' | awk '{print $1"\t"$9"\t"$4"\t"$5}' | sed 's/ID=//' >zj.gff`
#这里是把geneID加上物种标识，用于种间比对区分
#得先处理gff3文件，提取出自己想要的转录本，把基因iD加上染色体前缀
`grep -Ff longestid.txt dmy_gff3 > dmy_MCscanX.txt`
`awk '$3 == "mRNA" ' ../zj.gff3 | awk -F ";" '{print $1}' |awk '{print "zj_" $0}' |awk '{print $1"\t"$9"\t"$4"\t"$5"\t"$9}' |awk '{gsub(/ID=/,"zj("$1")"); print}' >zj_M.gff `
提取一个对照做mapping.txt,同时改好名字
`awk '{print $5"\t"$5}' zj_M.gff | sed 's/^zj([^)]*)//g' >mapping.txt`
这里就把基因的ID替换为mapping对应的处理好的有标识的基因id
`awk 'NR==FNR {split($1, id, "."); map[id[1]]=$2; next} /^>/ {split($1, header, "."); sub($1, ">"map[substr(header[1], 2)])} {print}' zj_mapping.txt zj_P450.fasta`
#使用全蛋白序列进行共线性分析
sed -e 's/Chr/zj/g' -e 's/scaffold/zj/g' zj.gff3 |sed 's/ID=/ID=zj_/g' >zj_filter.gff3
#染色体和geneID加上zj_前缀
sed 's/^>/>zj_/' zj.aa >zj_filter.aa
#注释行加上zj_前缀
#建库
cat sy_filter.aa zj_filter.aa >zy.aa
makeblastdb -in zy.aa -dbtype prot -out zy
#blast比对
nohup blastp -query zy.aa -db zy -out zy -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -num_threads 30 &
#这里是把geneID加上物种标识，用于种间比对区分
#得先处理gff3文件，提取出自己想要的转录本，把基因iD加上染色体前缀
`grep -Ff longestid.txt dmy_gff3 > dmy_MCscanX.txt`
`awk '$3 == "mRNA" ' ../zj.gff3 | awk -F ";" '{print $1}' |awk '{print "zj_" $0}' |awk '{print $1"\t"$9"\t"$4"\t"$5"\t"$9}' |awk '{gsub(/ID=/,"zj("$1")"); print}' >zj_M.gff `
提取一个对照做mapping.txt,同时改好名字
`awk '{print $5"\t"$5}' zj_M.gff | sed 's/^zj([^)]*)//g' >mapping.txt`
这里就把基因的ID替换为mapping对应的处理好的有标识的基因id
`awk 'NR==FNR {split($1, id, "."); map[id[1]]=$2; next} /^>/ {split($1, header, "."); sub($1, ">"map[substr(header[1], 2)])} {print}' zj_mapping.txt zj_P450.fasta`
#种间共线性全蛋白组
#先改蛋白序列的名字，在转录本g.t前面加物种标识
sed 's/g[0-9]/sy_&/g' sy_rmTE.aa >sy.aa
#gff3文件同理
后续处理方式同种内法分析#计算染色体长度
seqkit fx2tab -l -n -i Oc.fasta | awk '{print $1"\t"$2}' >Oc.len
#
awk '{print "chr\t-\t"$1"\t"$1"\t0\t"$2"\tchr"NR}' Oc.len > Oc.circos.txt