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

关键词

Python | C++ | 组学数据 | 阵列 | 序列 | 可视化 | 量化 | 基因 | 质量 | DNA | RNA | 遗传变异 | 蛋白质 | 氨基酸 | 血红蛋白 | 模式匹配 | 文本 | 正则表达 | 概率 | 比对 | 距离 | 密码子 | 频率 | 数据库 | 原子 | 细胞

Python和C++组学数据DNA和RNA | Notion亚图跨际 on Notion

🏈page指点迷津 | Brief

🎯要点

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;
 }

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