iso-seq bed 格式数据来源
同https://www.jianshu.com/p/5c8fe8c2d209
stop codon注释来源
Schizosaccharomyces_pombe.ASM294v2.45.gff3
stop codon bed6 file生成步骤
step1.输出wtf,tf transposon ID 列表，用于过滤wtf和 tf transposon
less Schizosaccharomyces_pombe.ASM294v2.45.gff3.gz|awk '{if($2=="PomBase" && $3~"gene") print $0;}'|grep wtf|gawk 'match($0, /ID=gene:(.*);Name=(.*);biotype=/, a) {print a[2]"\t"a[1];}' >wtf_PomBase_ID.txt
less Schizosaccharomyces_pombe.ASM294v2.45.gff3.gz|awk '{if($2=="PomBase" && $3~"gene") print $0;}'|grep Tf|gawk 'match($0, /ID=gene:(.*);Name=(.*);biotype=/, a) {print a[2]"\t"a[1];}' >Tf_PomBase_ID.txt
Note:wtf11-antisense-1,tfg,tfa等基因也会匹配上，从列表中删除。
step 2.抽取CDS坐标,并过滤wtf(注释不准确),Tf
less Schizosaccharomyces_pombe.ASM294v2.45.gff3.gz |awk '{if($2=="PomBase" && $3=="CDS") print $0;}'|gawk 'match($0, /:pep;Parent=transcript:(.*)\.[0-9]*;protein_id=/, a) {print $1"\t"$4"\t"$5"\t"$7"\t"a[1];}'|perl /data/XYH/script/removerow.pl wtf_PomBase_ID.txt 2 - 5 Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.txt
perl /data/XYH/script/removerow.pl Tf_PomBase_ID.txt 2 Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.txt 5 Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.txt
step 3.将同一个基因的不同CDS连接并输出头尾坐标
cat Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.txt | sort -k5,5 -k1,1 -k4,4 -k2,2n |perl CDS_link.pl - Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.link.txt
step 4.生成cds bed6 file
cat Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.link.txt|awk '{ print $1"\t"$2-1"\t"$3"\t"$5"\t.\t"$4;}' |sort -k1,1 -k2,2n > Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.sort.bed
Note:上述五个步骤路径：/data/XYH/snapgene_gbk_generator/work.sh
step 6. 生成stop codon bed6 file
cat /data/XYH/snapgene_gbk_generator/Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.sort.bed |awk '{if($6=="+")print $1"\t"$3-3"\t"$3"\t"$4"\t"$5"\t"$6;else print $1"\t"$2"\t"$2+3"\t"$4"\t"$5"\t"$6;}' >Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.stop_codon.bed
wtf stop codon注释来源
因为上述gtf中不包含wtf stop codons,因此直接使用我们自己的注释
文件路径:
/data/XYH/reference/wtf1-25_gene.bed
输出bed6格式：
cat /data/XYH/reference/wtf1-25_gene.bed |awk '{if($6=="+")print $1"\t"$3-3"\t"$3"\t"$4"\t"$5"\t"$6;else print $1"\t"$2"\t"$2+3"\t"$4"\t"$5"\t"$6;}' >wtf.stop_codon.bed
用bedtools计算跨过每个stop codon 的分子数
bedtools intersect -s -a Schizosaccharomyces_pombe.ASM294v2.45.gff3.PomBase_CDS.nowtf.noTf.stop_codon.bed -b meiosis_${n}h.isoseq.bed -wo |cut -f 4 |sort |uniq -c|awk '{print $2"\t"$1;}'>meiosis_${n}h.isoseq.overlap_with_stopcodon.count.txt
wtf 的计算同上
参数 -s 控制reads 和stop codon 在同一链上
用add_feature.pl将几个sample整理成一个表
echo "gene" >gene.list
cat Schizosaccharomyces_pombe.ASM294v2.47.stop_codon.bed|cut -f 4 >>gene.list
perl /data/XYH/script/add_feature.pl -i gene.list -c1 1 -p meiosis_\*h.isoseq.overlap_with_stopcodon.count.txt -c2 1 -c3 2 -title -s .isoseq.overlap_with_stopcodon.count.txt -o meiosis_isoseq.overlap_with_stopcodon.count.txt
增加基因名
less -S /data/XYH/reference/Schizosaccharomyces_pombe.ASM294v2.47.gtf.gz |awk '{if($9=="gene_id" && $11=="gene_name") print $0;}' |awk 'BEGIN{RF="\t";}{match($10,/"(.+)";/,id);match($12,/"(.+)";/,name);print id[1]"\t"name[1];}' >ID_and_name.txt
perl /data/XYH/script/add_feature.pl -i meiosis_isoseq.overlap_with_stopcodon.count.txt -c1 1 -p ID_and_name.txt -c2 1 -c3 2 -o meiosis_isoseq.overlap_with_stopcodon.count.addname.txt
