数据：
基因组文件： reference genome file(pombe.ASM294v2.29.dna.toplevel.fa), JB4.raven.fasta, JB864.raven.fasta(PacBio assemblies by SuoFang)
工具：
blast, RepeatModeler, RepeatMasker, bedtools, homemade perl/shell scripts
第一部分：用blastn搜索outerrepeat
第1步：用RepeatModeler 查找基因组上重复序列
#建库：
~/Software/RepeatModeler-2.0.1/BuildDatabase -name ref -dir ./genome
#Note: -name 库命名参数 ./genome 文件夹中存放基因组文件
#运行RepeatModeler
~/Software/RepeatModeler-2.0.1/RepeatModeler -database ref -pa 8 -LTRStruct
#Note:（from GuoSong）两个核心输出文件
1） <database_name>-families.fa : Consensus sequences ie. "ref-families.fa"
"consensus sequence"文件可以用作 RepeatMasker 的自定义library
2） <database_name>-families.stk : Seed alignments ie. "ref-families.stk"
"seed alignment"文件是Dfam数据库适配的"Stockholm"格式，可以直接上传Dfam数据库。
#重命名"ref-families.fa"中的序列名，方便后续分析 e.g."rnd-2_family-77#LTR/Gypsy ( Recon Family Size = 23, Final Multiple Alignment Size = 22 )" ——> “rnd-2_family-77”
cat ref-families.fa | awk '{if($0~/^>/){match($0,/>(.+)#/,a); print ">"a[1];}else{print $0;}}'|grep -v '^$' >ref-families.rename.fa
第2步：通过blast筛选ref-families.rename.fa中属于LTR和Tf的序列
#blast
makeblastdb -in /data/XYH/assemblies_by_wolf_et/LTR.fa -out LTR -dbtype nucl
blastn -query ref-families.rename.fa -db LTR -outfmt 7 -out ref-families.to_LTR.out
#blast 结果处理：将比对到LTR的序列删选出来并调整至方向一致
perl /data/XYH/script/blast_1st_hits.pl ref-families.to_LTR.out ref-families.to_LTR.1sthit.out
cat ref-families.to_LTR.1sthit.out|awk '{if($9<$10) print $1"\t+";else print $1"\t-";}'>ref-families.strand.txt
perl /data/XYH/script/fasta_length.pl ref-families.rename.fa |perl /data/XYH/script/add_feature.pl -i - -c1 1 -p ref-families.strand.txt -c2 1 -c3 2 -o ref-families.reverse.tmp
cat ref-families.reverse.tmp |awk '{if($3!="--")print $1"\t1\t"$2"\t"$3"\t"$1;}'>ref-families.reverse.txt
perl /data/XYH/script/03_extract_fasta_from_ref.pl ref-families.rename.fa ref-families.reverse.txt ref-families.reverse.fa
#加入标准LTR序列，生成alignment file
cat /data/XYH/assemblies_by_wolf_et/LTR.fa ref-families.reverse.fa |mafft --quiet - >ref.repeatmodeler-families.LTR.mfa
jalview打开，删除箭头所指序列
ref.repeatmodeler-families.LTR.mfa
Note:reference 中删除的序列比对到pomTelo2R,JB4.raven和JB864.raven删除的序列均有一段blast到wtf,是由于S.pombe的错误注释,将一段LTR包含在了其中
第3步：加入标准Tf CDS，LTR，生成RepeatMasker library
cat /data/XYH/reference/Tf_CDS.fa /data/XYH/assemblies_by_wolf_et/LTR.fa ref-families.reverse.fa |awk 'BEGIN{RS=">";}{if($1!="rnd-1_family-11" && $1!="") print ">"$0;}'|sed 's/rnd-.*_family/ref_queryLTR/' >ref.query_LTR.fa
#Note:将rnd-.*_family的序列名替换是为了区分来自不同基因组的query LTR,非必须操作
第4步：运行RepeatMasker，筛选比对到LTR或Tf的区域，生成bed6格式文件
RepeatMasker -lib ../ref.query_LTR.fa pombe.ASM294v2.29.dna.toplevel.fa
cat pombe.ASM294v2.29.dna.toplevel.fa.out | awk '{if($10~"ref" || $10~"Tf"){print $0;}}'|awk '{if($9=="+"){str="+";}else{str="-";} print $5"\t"$6-1"\t"$7"\trepeatmasker"NR"\t.\t"str;}' >ref.repeatmasker_LTR.bed
第5步：用bedtools maskfasta 将基因组中第4步比对到LTR的区域mask，并将连续N替换为单个N
cat reference_LTR_mask.fa |awk 'BEGIN{RS=">";FS="\n";}{if($1!=""){seq="";for (i=2; i<NF; i++){ seq=(seq""$i);} print ">"$1"\n"seq;}}'|sed 's/N\+/N/g' >reference_LTR_mask.single_letter.fa
#Note:将连续N替换为单个N是为了更容易找到可能被LTR隔开的outerrepeat
第6步：用blastn搜索outerrepeat,-task blastn，表示somewhat blast,其他均为默认参数
makeblastdb -in reference_LTR_mask.single_letter.fa -out ref_mask -dbtype nucl
for id in wtf19 wtf24 wtf4
do
blastn -query ../../outerrepeat_JB22_${id}-dw.fa -db ref_mask -task blastn -outfmt 7 -out ref_mask.query_${id}-dw.out
done
第7步：注释hits与wtf的位置关系
#在LTR masked,并将连续N替换为单个N的基因组中定位wtf
blastn -query REF.wtf.routine.fa.reverse.fa -db ref_mask -outfmt 7 -out ref_mask.wtf.out
#生成wtf bed6文件
perl /data/XYH/script/blast_perfect_hits.pl REF.wtf.routine.fa.reverse.fa ref_mask.wtf.out |awk '{if($9<$10) print $2"\t"$9-1"\t"$10"\t"$1"\t.\t+";else print $2"\t"$10-1"\t"$9"\t"$1"\t.\t-";}' |sort -k1,1 -k2,2n>ref_mask.wtf.bed
for id in wtf19 wtf24 wtf4
do
#生成outerrepeat hits bed6文件
cat ref_mask.query_${id}-dw.out |awk '{if($1!~/^#/) print $0;}' |awk '{if($9<$10) print $2"\t"$9-1"\t"$10"\t"$1"\t.\t+";else print $2"\t"$10-1"\t"$9"\t"$1"\t.\t-";}' |sort -k1,1 -k2,2n>ref_mask.query_${id}-dw.bed
#计算每个wtf与距离最近的hits的距离
bedtools closest -k 2 -s -D a -a ref_mask.wtf.bed -b ref_mask.query_${id}-dw.bed |awk '{if($13>=-1000 && $13<=1000) print $0;}'|awk '{if($13>0){posi="-down"}else{posi="-up"} split($4,a,"_"); print $7"\t"$8+1"\t"$9"\t"$12"\t"a[2]posi"("$13"bp)"}'|sort -k5,5 |sort -k1,1 -k2,2 -k3,3 -u >ref_mask.query_${id}-dw.distowtf.txt
#Note: -k 2 此参数非常重要，当基因上下游都存在outerrepeat时i，只有加这个参数才能全部注释
-s 参数要求hits与链同一方向，
-D a 参数使在wtf上游的hits与wtf的距离值为负数利用这一信息，可以添加注释-up/-down
#将hits相对于wtf的位置信息添加到blast结果中
perl /data/XYH/script/addanno_to_hits.pl -i ref_mask.query_${id}-dw.out -f ref_mask.query_${id}-dw.distowtf.txt -c1 1,2,3 -c2 5 -o ref_mask.query_${id}-dw.withanno.out
#抽取wtf上下游的outerrepeat hits
perl /data/XYH/script/03_extract_fasta_from_ref.pl reference_LTR_mask.single_letter.fa ref_mask.query_${id}-dw.distowtf.txt ref_mask.query_${id}-dw.outerrepeat.fa
#抽取序列过程生成的.fai文件最好每次都删除，因为这个文件无法被覆盖
rm reference_LTR_mask.single_letter.fa.fai
#生成outerrepeat alignment
mafft ref_mask.query_${id}-dw.outerrepeat.fa >ref_mask.query_${id}-dw.outerrepeat.mfa
done
最终结果
genome:JB864.raven.fasta，query:wtf19-dw
JB864.raven_mask.query_wtf19-dw.distowtf.mfa
blast result with position information
e.g.wts0403a-down(485bp)表示wts0403a下游outerrepeat在离stop codon 485bp处
第二部分：用repeatmasker搜索outerrepeat
从第一部分中LTR masked,并将连续N替换为单个N的基因组出发，这次将JB22_wtf4-down,JB22-wtf19-down,JB22_wtf24-down 3个queries放在一个文件中做为Repeatmasker library
第一步：运行Repeatmasker,筛选hits,生成bed6文件
#运行Repeatmasker
RepeatMasker -lib ../../3queries.fa reference_LTR_mask.single_letter.fa
#生成bed6文件
cat reference_LTR_mask.single_letter.fa.out |awk '{if($10~"JB22") print $0;}' |awk '{if($9=="+"){str="+";}else{str="-";} print $5"\t"$6-1"\t"$7"\trepeatmasker"NR"\t.\t"str;}' |sort -k1,1 -k2,2n|bedtools merge -s -c 6 -o distinct -i - |awk '{print $1"\t"$2"\t"$3"\touterrepeatbyRM"NR"\t.\t"$4}'>reference_RMouterrepeat.bed
第2步：注释hits与wtf的位置关系，与第一部分计算每个wtf与距离最近的hits的距离操作相同
bedtools closest -s -D a -a ref_mask.wtf.bed -b reference_RMouterrepeat.bed |awk '{if($13>=-1000 && $13<=1000) print $0;}'|awk '{if($13>0){posi="-down"}else{posi="-up"} split($4,a,"_"); print $7"\t"$8+1"\t"$9"\t"$12"\t"a[2]posi"("$13"bp)"}'|sort -k5,5 |sort -k1,1 -k2,2 -k3,3 -u >reference_RMouterrepeat.distowtf.txt
第3步：抽取outerrepeat序列并生成alignment
perl /data/XYH/script/03_extract_fasta_from_ref.pl reference_LTR_mask.single_letter.fa reference_RMouterrepeat.distowtf.txt reference_RMouterrepeat.distowtf.fa
mafft reference_RMouterrepeat.distowtf.fa >reference_RMouterrepeat.distowtf.mfa
