工具：

bedtools, a perl script(abundant_tss_count.pl)

数据：

1."/data/XYH/reference/wtf1-25_gene.bed":  bed6 file recording positions of  25 wtf genes in the reference genome

2."meiosis_*h.isoseq.bed":  transformated from isoseq "rep2"  bam files with bedtools bamtobed

3.“up_conserved_up.type.txt” ：copied from the Email  2020年5月19日 (周二) 11:40 "Subject: upstream sequence types and Iso-Seq" from Dr.Du

流程：

step1: 查找落在每个wtf gene上的iso-seq reads,并计算read 第一个碱基相对于wtf conserved_up第一个碱基的位置

bedtools intersect -s -f 0.5 -a /data/XYH/reference/wtf1-25_gene.bed -b meiosis_${n}h.isoseq.bed -wo |awk '{if($6 =="+") {print $4"\t"$8-$2+1;} else {print $4"\t"$3-$9+1;}}' >meiosis_${n}h.wtf.tss.txt0

重要参数：

-s  控制reads与基因同一链方向

-f 0.5 read 与基因overlap的长度占基因总长度（从conserved_up到stop codon）50%以上

awk这部分命令用于计算每条read的转录起始位置

step1 result:

step2:用脚本abundant_tss_count.pl统计每个wtf 基因的reads数目，查找有多条reads一致的TSS sites，并计算在这一位置起始的reads占这个基因total reads的百分比，并输出上游序列类型信息以及以上的统计信息

重要参数设置

1.对于total reads <10条的基因不统计

2.对于某个位置,如果在该位置起始的reads 数占该基因total reads 10%以上或者表达量较高时，虽然在该位置起始的reads数目不到10%，但是 绝对数目>10时，报告这个位置，作为一个tss site,并输出在该位置起始的reads数目占total reads的百分比信息

step2 result