今天运行一个程序的时候，出现了一个让我惊讶的报错，
报错
照理说，基因组上的序列应该只有ATCG，以及用来填充gap的N才对，为啥会出现Y呢？于是我用grep又去搜索了下，出现了更多让我看不懂的字符。
离奇现象
为了一探究竟，我展开了一番搜索。由于这个现象不知道怎么描述，于是就随便试试，当然都没有成功。突然脑子闪过一个想法，这不会碱基命名中的用来表示哪些不好判断的序列呀，因为我想到在对基因型分型的时候，如果只能确定这个基因型不是AA，但是有可能是AB,BB，那么就会用另一个字母进行表示。
于是我想到了IUPAC命名法，最后我找到了下表
| IUPAC nucleotide code | Base |
|---|---|
| A | Adenine |
| C | Cytosine |
| G | Guanine |
| T (or U) | Thymine (or Uracil) |
| R | A or G |
| Y | C or T |
| S | G or C |
| W | A or T |
| K | G or T |
| M | A or C |
| B | C or G or T |
| D | A or G or T |
| H | A or C or T |
| V | A or C or G |
| N | any base |
| . or - | gap |
解决了心中的一个疑问。那么下一个问题就是如何处理这些非ATCGN的字符呢？我当然使用最简单粗暴的方法，就是把RYSWKMBDHV全部替换成N
tr RYSWKMBDHVryswkmbdhv N < /data/reference/genome/TAIR10/Athaliana.fa | sed -e 's/>CNN/>Chr/' -e 's/>ChrN/>ChrM/' -e 's/>ChrCN/>ChrCh/' > r
ef/Athaliana.fa
# 由于把Chr也替换成CNN, 因此要替换回来
网友评论