Python和C++组学数据DNA和RNA

🎯要点

1. Python和C++计算：
   1. 🎯生物计算：🖊分析和编辑序列文件 | 🖊微阵列数据分析：阵列图像转换为以二进制或文本格式保存的显着强度或可量化值 、可视化数据分布质量控制、分位数归一化比较数据、差异基因表达分析 | 🖊遗传或物理共表达基因的聚类图 | 🖊基因富集分析生物过程和分子途径 | 🖊单核苷酸多态性分析：药物反应与遗传变异的关联分析。🎯处理方式：🎯解析：🖊质谱数据整合到代谢途径 | 🖊核苷酸序列或氨基酸（蛋白质）序列文件解析 | 🖊序列数据库解析 | 🎯搜索：🖊RNA 序列转为蛋白质序列 | 🖊核苷酸序列或氨基酸（蛋白质）序列文件填入字典，定义对应的键和值 | 🖊创建蛋白质序列循环预测器 | 🖊从蛋白质数据库文提取氨基酸序列 | 🎯分类：🖊序列同一性百分比对输出排序 | 🖊对血红蛋白条目排序 | 🎯模式匹配和文本挖掘：🖊蛋白质正则表达式转换为 Python 正则表达式 | 🖊寻找基因组序列转录因子结合位点 | 🎯代码计算：🖊从核苷酸序列或氨基酸（蛋白质）序列文件文件加载蛋白质序列，计算每列的间隙分数 | 🖊创建与给定序列具有相同氨基酸或核苷酸组成的随机序列 | 🖊创建具有相同概率的核苷酸随机序列 | 🖊解析多个序列比对 | 🖊从多序列比对计算共识序列 | 🖊计算系统发育树节点之间的距离 | 🖊核苷酸序列中的密码子频率 | 🖊解析维也纳格式的 RNA 二维结构 | 🖊解析生物数据库 XML 输出和系统生物学标记语言文件 | 🖊执行生物数据库 | 🖊将蛋白质数据库文件拆分为蛋白质数据库链文件 | 🖊查找蛋白质数据库结构中两个最接近的原子 | 🖊提取两个蛋白质链之间的接口 | 🖊构建同源模型 | 🖊RNA 三维同源建模 | 🖊从 3D 结构计算 RNA 碱基对 | 🖊丝氨酸蛋白酶催化三联体。
2. 分析工具：🎯分层交错布隆结构，快速查找生物序列 | 🎯量化大型测序数据集 | 🎯距离编辑精确成对序列比对 | 🎯高效的细胞计数分析

🍇C++和脚本预处理生物序列

预处理的一个主要用例是在 FASTA 和 FASTQ 之间转换数据。

awk '{if(NR%4==1){print ">" substr($1, 2)}else if(NR%4==2){print}}' FILE.fastq > FILE.fasta
awk '{if(NR%4==1){print ">" substr($1, 2)}else if(NR%4==2){print}}' FILE.fastq | fold > FILE.fasta
awk '{if(NR==1 && $0 ~ /^>/){print}else if (NR!=1 && $0 ~ /^>/){print "\n" $0}else {printf $0}} END {printf "\n"}' FILE.fasta > NEW_FILE.fasta

重命名 FASTQ 读取

awk '/^@EXISTING_NAME_REGEX/{print "@TARGET_NAME_" ++i; next}{print}' FILE.fastq > NEW_FILE.fastqv

读取GZ文件

zcat FILE.tar.gz | head -nNO_READS > FILE.fastq

为了绘制读数长度的直方图和其他实验（例如过滤短读数），我们可以执行以下操作。请注意，我们需要对长读应用程序进行此类预处理。

awk "{{ if(NR%4==2) print length($1) }}" READS.fastq> LENGTHS.txt

将双端读取分成两个文件

awk 'c-->0;$0~s{if(b)for(c=b+1;c>1;c--)print r[(NR-c+1)%b];print;c=a}b{r[NR%b]=$0}' b=0 a=3 s="/1" READS.fastq > R1.fastq
awk 'c-->0;$0~s{if(b)for(c=b+1;c>1;c--)print r[(NR-c+1)%b];print;c=a}b{r[NR%b]=$0}' b=0 a=3 s="/2" READS.fastq > R2.fastq

在此计算中，我们以二进制形式对 k-mers 进行编码。 A-00、C-01、T-10 G-11 使用以下方法。

CHAR >> 1 & 3;

我们在下面的函数中使用上述方法。 这里我们传递一个字符串引用和间隔来计算所需长度的 k-mers。 然而，由于我们使用 u_int64_t，因此我们仅限于 32 聚体。 对于除此之外的任何内容，请尝试使用位集或布尔向量。

u_int64_t substr_kmer_int(string &str, int start, int length)
{
    u_int64_t val = 0;

    for (size_t i = start; i < start + length; i++)
    {
        val = val << 2;
        val += (str[i] >> 1 & 3);
    }

    return val;
}

然后我们使用下面的一组位运算来获得反向补码。

u_int64_t revComp(u_int64_t x, size_t sizeKmer)
{
    u_int64_t res = x;

    res = ((res >> 2 & 0x3333333333333333) | (res & 0x3333333333333333) << 2);
    res = ((res >> 4 & 0x0F0F0F0F0F0F0F0F) | (res & 0x0F0F0F0F0F0F0F0F) << 4);
    res = ((res >> 8 & 0x00FF00FF00FF00FF) | (res & 0x00FF00FF00FF00FF) << 8);
    res = ((res >> 16 & 0x0000FFFF0000FFFF) | (res & 0x0000FFFF0000FFFF) << 16);
    res = ((res >> 32 & 0x00000000FFFFFFFF) | (res & 0x00000000FFFFFFFF) << 32);
    res = res ^ 0xAAAAAAAAAAAAAAAA;

    return (res >> (2 * (32 - sizeKmer)));
}

计算 K-mer 直方图

我们可以扩展上述函数来获取序列的 K-mer 直方图，这在生物信息学分析中非常常用。在第一个函数中，我们获得给定 k 值的所有可能的 K-mer。这是为了通过仅更新查找表来使计算更快。

unordered_map<int, double> kSizeMers(int size)
{
    unordered_map<int, double> kmers;

    for (u_int64_t i = 0; i < (int)pow(4, size); i++)
    {
        kmers[min(i, revComp(i, size))] = 0;
    }

    return kmers;
}

vector<double> getKmersProfile(string &line, int size, bool normalize, unordered_map<int, double> &KMERS)
{
    double total = 0;
    long len = 0;
    u_int64_t val = 0;
    unordered_map<int, double> kmers;
    kmers.insert(KMERS.begin(), KMERS.end());
    vector<double> stats(kmers.size(), 0);

    for (int i = 0; i < (int)line.length(); i++)
    {
        if (!(line[i] == 'A' || line[i] == 'C' || line[i] == 'G' || line[i] == 'T'))
        {
            val = 0;
            len = 0;
            continue;
        }

        val = (val << 2);
        val = val & (int)(pow(2, 2 * size) - 1);
        val += (line[i] >> 1 & 3);
        len++;

        if (len == size)
        {
            // use val as the kmer for counting
            len--;
            kmers[min(val, revComp(val, size))]++;
            total++;
        }
    }

    int i = 0;
    for (unordered_map<int, double>::iterator it = kmers.begin(); it != kmers.end(); ++it)
    {
        if (normalize)
        {
            stats[i] = it->second / max(1.0, total);
        }
        else 
        {
            stats[i] = it->second;
        }
        i++;
    }

    return stats;
}

在第二个函数中，我们迭代整个字符串并在查找表的副本中记录最低补码。我们进行复制，以便可以在多线程场景中重复使用相同的原始查找表。

参阅一：计算思维

参阅二：亚图跨际