get mutated protein sequence from all mutation(snp and indels)
requirements
pandas
vcf
loguru
tqdm
biopython
bcftools
seqkit
csvkit
sh ./download.sh 109 # download seqnence GCRh38 109 release from ensembl
1,109版本中有121766个转录本其中99812个转录本是以三个常见的终止密码子停止转录的,但是依然有21954个密码子不是以终止密码子停止.
在这21954个转录本中,切取最后三个碱基统计:
GGC | AGT | ACG | TAC | GGA | TCC | GAC | ACC | CAC | AAC | GCC | AGC | GTC | TGG | AGA | ATA | ATC | GGT | CTC | AGG | ATG | TTG | TTC | TGT | GCT | ACA | GTA | TGC | GCG | GGG | GCA | ACT | TTA | TCG | GAG | GTG | CAG | CTA | TCT | AAA | CCT | TAT | GAT | TCA | AAT | CGG | ATT | CAA | CCC | CCA | CTG | CAT | CGT | CGA | CTT | TTT | AAG | CCG | GAA | CGC | GTT |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
276 | 208 | 148 | 209 | 512 | 368 | 257 | 361 | 268 | 258 | 371 | 339 | 206 | 652 | 613 | 168 | 315 | 273 | 301 | 529 | 494 | 315 | 254 | 273 | 381 | 347 | 125 | 312 | 157 | 425 | 394 | 275 | 131 | 136 | 754 | 423 | 1205 | 156 | 353 | 662 | 590 | 165 | 254 | 351 | 260 | 266 | 223 | 420 | 524 | 487 | 771 | 320 | 119 | 143 | 393 | 310 | 1031 | 275 | 529 | 161 | 158 |
同样,在终止位点向后移动三个碱基作为统计:
GGC | AGT | ACG | TAC | TGT | TGG | AGC | GTC | AGA | ATA | ATC | GGT | CTC | AAC | ACC | AGG | ATG | TTG | TTC | TGA | GCT | GGA | ACA | GTA | TGC | GCG | GGG | GCA | ACT | TTA | TCG | GAG | GAC | CTA | TCC | TCT | CAG | TTT | GTG | CCT | GAA | TCA | AAT | CGG | TAT | AAA | ATT | CAA | CCC | CCA | CTG | CAT | CAC | CGT | CGA | CTT | AAG | CGC | TAG | GAT | CCG | GCC | GTT | TAA | NNN | NNG |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
318 | 212 | 175 | 221 | 341 | 655 | 396 | 242 | 432 | 153 | 309 | 276 | 389 | 250 | 1508 | 453 | 362 | 377 | 361 | 201 | 361 | 610 | 376 | 139 | 372 | 138 | 360 | 375 | 431 | 186 | 121 | 477 | 210 | 166 | 381 | 495 | 841 | 272 | 327 | 642 | 377 | 333 | 214 | 222 | 147 | 404 | 229 | 320 | 381 | 437 | 649 | 278 | 359 | 124 | 131 | 434 | 621 | 154 | 102 | 192 | 229 | 385 | 214 | 86 | 16 | 5 |
结论:大部分转录本(82%)是以终止子停止转录蛋白的,但是也有部分(18%)的转录本不以标准的密码子就结束了转录,在处理时,我们对于以标准密码子结束的转录本对翻译结束位点延长一些长度(如5000bp),但是对于不以标准的密码子结束的转录本我们不进行延长。
# 替换非标准氨基酸
sed '/^[^>]/s/[BJOUZ]/X/g' all_mutated.fa > all_mutated_filter.fa
#计算embedding
alias esmExt="python /data/wenyuhao/code/someProject/esm/scripts/extract.py esm2_t36_3B_UR50D"
muteFa=/data/wenyuhao/55/ensembl/ensembl98/MS/MS2/result/all_mutated_filter.fa
esmExt $muteFa all_mutated_filter --repr_layers 0 35 36 --include mean per_tok