Ribo-seq分析必看文献 | 知识（五）：RiboTaper
2019-03-07 本文已影响3人
热衷组培的二货潜
一、echo "Extracting gene names + biotypes from gtf..."
一、得到基因的对应的注释信息，基因ID|基因类型|基因名称
#############################################################################################
##1.1、得到 `基因ID|基因类型|基因名称 ` 的文件
#############################################################################################
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") ## 如果以gene_id开头的，则记录当前行号x
for (y=1;y<=NF;y++) if ($y~/gene_type|gene_biotype/) ## 以gene_type或者gene_biotype开头的行号
for (z=1;z<=NF;z++) if ($z~"gene_name") ## 以gene_name开头的行号
print $(x+1) "\t" $(y+1) "\t" $(z+1)}' RiboTaper_test.gtf | sort | uniq | sed 's/;//g' | sed 's/"//g' > gene_name_type
`
ENSG00000008128.18 protein_coding CDK11A
ENSG00000008130.11 protein_coding NADK
ENSG00000067606.11 protein_coding PRKCZ
ENSG00000078369.13 protein_coding GNB1
ENSG00000078808.12 protein_coding SDF4
ENSG00000107404.13 protein_coding DVL1
ENSG00000116151.9 protein_coding MORN1
ENSG00000127054.14 protein_coding CPSF3L
ENSG00000130762.10 protein_coding ARHGEF16
ENSG00000131584.14 protein_coding ACAP3
`
#############################################################################################
##1.2、 从gtf文件中得到没有指定type的行
#############################################################################################
less gene_name_type | cut -f 1 | grep -Fvf - RiboTaper_test.gtf |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (z=1;z<=NF;z++) if ($z~"gene_name")
print $(x+1) "\t" "no_biotype" "\t" $(z+1)}' | sort | uniq | sed 's/;//g' | sed 's/"//g' > gene_name_notype
# http://man.linuxde.net
grep :
-F 将范本样式视为固定字符串的列表。
-f<范本文件> 指定范本文件，其内容有一个或多个范本样式，让grep查找符合范本条件的文件内容，格式为每一列的范本样式。
-v 反转查找。
#############################################################################################
##1.3、 从gtf文件中得到没有gene名的行
#############################################################################################
less gene_name_type | cut -f 1 | grep -Fvf - RiboTaper_test.gtf |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~/gene_type|gene_biotype/)
print $(x+1) "\t" $(y+1) "\t" "no_name"}' | sort | uniq | sed 's/;//g' | sed 's/"//g' > gene_noname_type
#############################################################################################
##1.4、 合并gene_name_type gene_name_notype gene_noname_type三个文件到一个文件中
#############################################################################################
cat gene_name_type gene_name_notype gene_noname_type > gene_annot_name_pre
#############################################################################################
##1.5、 从gtf文件中提取既没有指定基因name的又没有指定基因type的行
#############################################################################################
less gene_annot_name_pre | cut -f 1 | grep -Fvf - RiboTaper_test.gtf |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
print $(x+1) "\t" "no_biotype" "\t" "no_name"}' | sort | uniq | sed 's/;//g' | sed 's/"//g' > gene_noname_notype
#############################################################################################
##1.5、 得到最终基因的注释文件
#############################################################################################
cat gene_annot_name_pre gene_noname_notype | sort | uniq > gene_annot_names
#############################################################################################
##1.6、 删除中间文件
#############################################################################################
rm gene_name_type gene_name_notype gene_noname_type gene_annot_name_pre gene_noname_notype
二、 echo "creating bed_files..."
#############################################################################################
##2.1、#TAKE CDS OF CCDS REGIONS
##2.1、得到CDS进一步的CCDS区域
#############################################################################################
if [ "$3" = true ] ; then
awk '{if($3=="CDS") print $0}' RiboTaper_test.gtf | grep ccdsid |
## 先筛选第三列为CDS的行，然后再筛选有ccdsid的行
awk '{ for (x=1;x<=NF;x++) if ($x~"^gene_id") print $1 "\t" $4-1 "\t" $5 "\t" "CCDS" "\t" $(x+1) "\t" $7 }' |
## 由于bed文件是从0开始的，所以在gtf文件的基础上要减去1，从而得到CCDS的坐标bed文件
sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > unique_ccds.bed
## 进行排序，将`;`和 `"`替换掉
less RiboTaper_test.gtf | grep ccdsid |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~"^transcript_id")
if($3=="exon") print $1 "\t" $4-1 "\t" $5 "\t" $(y+1)"\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_ccds_ccdsid.bed
fi
if [ "$3" = false ] ; then
## 没有ccdsid的注释gtf文件时候采用此命令
awk '{if($3=="CDS") print $0}' RiboTaper_test.gtf |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
print $1 "\t" $4-1 "\t" $5 "\t" "CCDS" "\t" $(x+1) "\t" $7 }' | sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > unique_ccds.bed
fi
#############################################################################################
##2.2、#STORE CCDS GENES
##2.2、得到CCDS相关基因的转录本名ID
#############################################################################################
less unique_ccds.bed | cut -f 5 | sort | uniq > genes_ccds
ENSG00000008128.18
ENSG00000008130.11
#############################################################################################
##2.3、#STORE COORDINATES CDS CCDS
##2.3、提取ccds文件的bed文件将原来的1到3列的信息”chr1 69090 70005“转变为”chr1_69090_70005“
#############################################################################################
less unique_ccds.bed | cut -f 1-3 | tr '\t' '_' > coords_ccds
chr1_69090_70005
chr1_367658_368594
#############################################################################################
##2.4、#TAKE ALL EXONS OF CCDS GENES
##2.4、不懂，看文章才清楚！！ 暂时放着
#############################################################################################
grep -Ff genes_ccds RiboTaper_test.gtf | awk '{if($3=="exon") print $0}' |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") print $1 "\t" $4-1 "\t" $5 "\t" "EXONCCDS" "\t" $(x+1) "\t" $7 }' |
sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > unique_exons_ccds.bed
#############################################################################################
##2.5、#STORE COORDINATES EXONS CCDS
##2.5、同2.3，将`xx xx xx`转换成`xx_xx_xx`的格式
#############################################################################################
less unique_exons_ccds.bed | awk '{print $1"_"$2"_"$3 "\t" $0}' > coords_unique_exons_ccds.bed
#############################################################################################
##2.6、#TAKE OUT CDS CCDS FROM EXONS CCDS
##2.6、挑选unique exon的 ccds bed坐标
#############################################################################################
grep -Fvf coords_ccds coords_unique_exons_ccds.bed | awk '{print $2 "\t" $3 "\t" $4 "\t" "EXONCCDS" "\t" $6 "\t" $7}' > unique_exons_ccds.bed
#############################################################################################
##2.7、#REMOVE COORDS
##2.7、筛除中间文件
#############################################################################################
rm coords_ccds coords_unique_exons_ccds.bed
三、 提取没有CCDS的基因的exon
#############################################################################################
## TAKE EXONS OF NONCCDS GENES
## 提取没有CCDS的基因的exon
#############################################################################################
grep -Fvf genes_ccds RiboTaper_test.gtf |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") if($3=="exon")
print $1 "\t" $4-1 "\t" $5 "\t" "EXONnonCCDS" "\t" $(x+1) "\t" $7 }' |
sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > unique_nonccds.bed
四、 提取没有CCDS的基因的exon
#############################################################################################
## #TAKE SEQUENCES, STRANDED INFO
## 提取序列
#############################################################################################
# bedtools中的说明书
# https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html
-tab Report extract sequences in a tab-delimited format instead of in FASTA format.
-fo Specify an output file name. By default, output goes to stdout.
-fi <input FASTA>
-bed <BED/GFF/VCF>
-name Use the “name” column in the BED file for the FASTA headers in the output FASTA file.
fastaFromBed -s -fi $genome_full -bed unique_ccds.bed -fo unique_ccds_seq.fa
# $genome_full 为已建立索引的基因组fasta序列
awk '{$4=$1"_"$2"_"$3"_"$4"_"$5"_"$6; print $0}' OFS="\t" unique_ccds.bed |
fastaFromBed -fi $genome_full -name -bed - -tab -fo unique_ccds_seq_name_tab
paste <(cut -f 1 unique_ccds_seq_name_tab| tr '_' '\t') <(cut -f 2 unique_ccds_seq_name_tab | sed 's/[A-Z]/& /g') > sequences_ccds
awk '{$4=$1"_"$2"_"$3"_"$4"_"$5"_"$6; print $0}' OFS="\t" unique_exons_ccds.bed |
fastaFromBed -fi $genome_full -name -bed - -tab -fo unique_exons_ccds_seq_name_tab
paste <(cut -f 1 unique_exons_ccds_seq_name_tab | tr '_' '\t') <(cut -f 2 unique_exons_ccds_seq_name_tab | sed 's/[A-Z]/& /g') > sequences_exonsccds
awk '{$4=$1"_"$2"_"$3"_"$4"_"$5"_"$6; print $0}' OFS="\t" unique_nonccds.bed |
fastaFromBed -fi $genome_full -name -bed - -tab -fo unique_nonccds_seq_name_tab
paste <(cut -f 1 unique_nonccds_seq_name_tab | tr '_' '\t') <(cut -f 2 unique_nonccds_seq_name_tab | sed 's/[A-Z]/& /g') > sequences_nonccds
fastaFromBed -s -fi $genome_full -bed unique_exons_ccds.bed -fo unique_exons_ccds_seq.fa
fastaFromBed -s -fi $genome_full -bed unique_nonccds.bed -fo unique_exons_nonccds_seq.fa
五、 提取找到的ORF的序列
#############################################################################################
## #CAT SEQUENCES TOGETHER FOR ORF FINDING
## 提取序列
#############################################################################################
cat unique_ccds_seq.fa unique_exons_ccds_seq.fa > unique_ccds_exonccds_seq.fa
六、 提取所有CDS的区域
#############################################################################################
## #make all CDS regions
## 提取CDS区域
#############################################################################################
less $genc_full | awk '{if($3=="CDS") print $0}' |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") print $1 "\t" $4-1 "\t" $5 "\t" "cds" "\t" $(x+1) "\t" $7 }' |
sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > all_cds.bed
七、 转录本信息
#############################################################################################
## #assembling transcript information...
## 组装转录本信息
#############################################################################################
#TAKE TRANSCR CCDS
grep -Ff genes_ccds $genc_full | awk '{ for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~"^transcript_id") if($3=="exon")
print $1 "\t" $4-1 "\t" $5 "\t" $(y+1) "\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_ccds.bed
八、#TAKE TRANSCR APPRIS CCDS
#############################################################################################
## #TAKE TRANSCR APPRIS CCDS
##
#############################################################################################
# if [ "$4" = true ] ; then
# use_appris = "true" or "false"
grep -Ff genes_ccds $genc_full | grep appris_prin |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") for (y=1;y<=NF;y++)
if ($y~"^transcript_id") if($3=="exon") print $1 "\t" $4-1 "\t" $5 "\t" $(y+1)"\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_ccds_appris_prin.bed
cut -f 5 transcr_exons_ccds_appris_prin.bed | sort | uniq > genes_appris_prin
grep -Ff genes_ccds $genc_full | grep -Fvf genes_appris_prin - |
grep appris | awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~"^transcript_id") if($3=="exon") print $1 "\t" $4-1 "\t" $5 "\t" $(y+1)"\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_ccds_appris_noprin.bed
cut -f 5 transcr_exons_ccds_appris_noprin.bed | sort | uniq > genes_appris_noprin
grep -Ff genes_ccds $genc_full | grep -Fvf genes_appris_prin - | grep -Fvf genes_appris_noprin |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~"^transcript_id") if($3=="exon") print $1 "\t" $4-1 "\t" $5 "\t" $(y+1)"\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_ccds_noappris_noprin.bed
cut -f 5 transcr_exons_ccds_noappris_noprin.bed | sort | uniq > genes_noappris_noprin
cat transcr_exons_ccds_appris_prin.bed transcr_exons_ccds_appris_noprin.bed transcr_exons_ccds_noappris_noprin.bed > transcr_exons_ccds_appris.bed
cat genes_appris_prin genes_appris_noprin genes_noappris_noprin > genes_ccds_appris
八、#TAKE TRANSCR NONCCDS
#############################################################################################
## #TAKE TRANSCR NONCCDS
##
#############################################################################################
grep -Fvf genes_ccds $genc_full | awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
for (y=1;y<=NF;y++) if ($y~"^transcript_id") if($3=="exon")
print $1 "\t" $4-1 "\t" $5 "\t" $(y+1)"\t" $(x+1) "\t" $7}'|
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k2,2n | uniq > transcr_exons_nonccds.bed
#start_stop_cds
less RiboTaper_test.gtf | awk '{for (x=1;x<=NF;x++)
if ($x~"^gene_id") if($3=="start_codon" || $3=="stop_codon")
print $1 "\t"$4-1 "\t"$5 "\t" $3 "\t" $(x+1) "\t"$7}' |
sed 's/;//g' | sed 's/"//g' | sort -k1,1 -k2,2g | uniq |
awk 'p{print $0 "\t" $2-p}{p=$2}' | tac | awk 'p{print $0 "\t" $2-p}{p=$2}' |
tac | awk '{if($NF<-100 || $(NF-1)>100) print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6}' > start_stops_FAR.bed
# tac 命令用于将文件已行为单位的反序输出，即第一行最后显示，最后一行先显示。
#make cds transcript coords
echo "Creating transcript cds coordinates from gtf..."
awk '{ for (x=1;x<=NF;x++) if ($x~"^transcript_id") if ( $3=="exon" || $3=="CDS" )
print $1 "\t" $3 "\t" $4 "\t" $5 "\t" $7 "\t" $(x+1)}' RiboTaper_test.gtf |
sed 's/"//g' | sed 's/;//g' | sort -k1,1 -k3,3g > exons_cds_all
$scripts_dir_full"/gtf_to_start_stop_tr.R"
··················
###script for creating transcript-level coordinates of CDS positions (start and stop) from a .gtf file
print(paste("--- extracting transcript-level CDS cordinates","---",date(),sep=" "))
exons_cds_all<-read.table("exons_cds_all",stringsAsFactors=F,header=F)
colnames(exons_cds_all)<-c("chr","type","start","end","strand","transcript_id")
tr_cds<-unique(exons_cds_all[exons_cds_all[,"type"]=="CDS","transcript_id"])
exons_cds_all2<-exons_cds_all[exons_cds_all[,"transcript_id"]%in%tr_cds,]
exons_cds_all2$length<-1+(exons_cds_all2$end-exons_cds_all2$start)
list_exons_cds_tr<-split.data.frame(x=exons_cds_all2,f=exons_cds_all2$transcript_id,drop=T)
list_coords<-list()
for(i in 1:length(list_exons_cds_tr)){
transcr<-tr_cds[i]
trascr_data<-list_exons_cds_tr[[transcr]]
strand<-trascr_data$strand[1]
exons_in_transcr<-trascr_data[trascr_data[,"type"]=="exon",]
if(strand=="-"){exons_in_transcr<-exons_in_transcr[dim(exons_in_transcr)[1]:1,]}
cds_in_transcr<-trascr_data[trascr_data[,"type"]=="CDS",]
if(strand=="-"){cds_in_transcr<-cds_in_transcr[dim(cds_in_transcr)[1]:1,]}
cumsumexons<-cumsum(exons_in_transcr$length)
revcumsumexons<-cumsum(rev(exons_in_transcr$length))
cumsumcds<-cumsum(cds_in_transcr$length)
st_cod<-cds_in_transcr[1,"start"]
if(strand=="-"){st_cod<-cds_in_transcr[1,"end"]}
end_cod<-(cds_in_transcr[dim(cds_in_transcr)[1],"end"])
if(strand=="-"){end_cod<-(cds_in_transcr[dim(cds_in_transcr)[1],"start"])}
st_ex<-which((st_cod>=exons_in_transcr$start & st_cod<=exons_in_transcr$end))
end_ex<-which((end_cod>=exons_in_transcr$start & end_cod<=exons_in_transcr$end))
nt_dist_start<-st_cod-exons_in_transcr[st_ex,"start"]
if(strand=="-"){nt_dist_start<-exons_in_transcr[st_ex,"end"]-st_cod}
if(st_ex>1){nt_dist_start<-nt_dist_start+cumsumexons[st_ex-1]}
nt_dist_stop<-exons_in_transcr[end_ex,"end"]-end_cod
if(strand=="-"){nt_dist_stop<-end_cod-exons_in_transcr[end_ex,"start"]}
if(end_ex<dim(exons_in_transcr)[1]){nt_dist_stop<-nt_dist_stop+revcumsumexons[dim(exons_in_transcr)[1]-end_ex]}
tr_len<-sum(exons_in_transcr$length)
start_coord<-nt_dist_start+1
stop_coord<-(tr_len-nt_dist_stop)+3
if(nt_dist_stop==0){stop_coord<-tr_len}
x<-data.frame(transcript_id<-transcr,start_tr<-start_coord,stop_tr<-stop_coord)
list_coords[[i]]<-x
}
coords<-do.call(what=rbind.data.frame,args=list_coords)
colnames(coords)<-c("transcript_id","start_tx","stop_tx")
write.table(coords,file="cds_coords_transcripts",row.names=F,col.names=F,quote=F,sep="\t")
print(paste("--- extracting transcript-level CDS cordinates, Done!","---",date(),sep=" "))
··················
#make cds frames
less RiboTaper_test.gtf | awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id") if($3=="CDS")
print $1 "_" $4-1 "_" $5 "_" "CCDS" "_" $(x+1) "\t" $8 "\t"$7 "\t" $5-($4-1)}' |
sed 's/;//g' | sed 's/"//g' | sort -k1,1 -k2,2 | uniq > frames_ccds
#take all exonic regions
less RiboTaper_test.gtf | awk '{if($3=="exon") print $0}' |
awk '{for (x=1;x<=NF;x++) if ($x~"^gene_id")
print $1 "\t" $4-1 "\t" $5 "\t" "exon" "\t" $(x+1) "\t" $7 }' |
sort -k1,1 -k2,2n | uniq | sed 's/;//g' | sed 's/"//g' > all_exons.bed
$scripts_dir_full"/genes_coor.R"
...............................................................................
all_ex<-read.table("all_exons.bed",stringsAsFactors=F,header=F)
spli<-split.data.frame(all_ex,f=all_ex$V5)
spli2<-lapply(spli,FUN=function(x){
minc<-min(x[,2])
maxc<-max(x[,3])
data.frame(chr=x[1,1],start=minc,end=maxc,le=maxc-minc,gene_id=x[1,5],strand=x[1,6],stringsAsFactors=F)
})
spli3<-do.call(what=rbind.data.frame,args=spli2)
write.table(file="genes_start_end",x=spli3,col.names=F,row.names=F,quote=F,sep="\t")
system("sort -k1,1 -k2,2n genes_start_end > genes_start_end.bed ")
system("rm genes_start_end")
..................................................................................
