在 NGS 数据中查找 CDR

我有数百万个fasta格式的序列，想要提取CDR（CDR1、CDR2和CDR3）。我只选择了一个序列作为示例，并尝试提取CDR1，但无法提取CDR1。

顺序：-'FYSHSAVTLDESGGGLQTPGGGLSLVCKASGFTFSSYGMMWVRQAPGKGLEYVAGIRNDA GDKRYGSAVQGRATISRDNGQSTVRLQLNNLRAEDTGTYFCAKESGCYWDSTHCIDAWGH GTEVIVSTGG'。

cdr1 开始于：-“VCKASGFTFS”，最多可更换 3 名，但第二名 C 是必须的。cdr1 结束于：-'WVRQAP'，最多可替换两名，但第三名的 R 是必须的。

提取的 cdr1 应该是SYGMM

def cdr1_in(cdr_in): #VCKASGFTFS
    pin=0
    max_pin=3       
    
    if cdr[1]!='C':
        pin+=1
    if cdr[0]!='V':
        pin+=1
    if cdr[2]!='K':
        pin+=1
    if cdr[3]!='A':
        pin+=1    
    if cdr[4]!='S':
        pin+=1
    if cdr[5]!='G':
        pin+=1
    if cdr[6]!='F':
        pin+=1
    if cdr[7]!='T':
        pin+=1    
    if cdr[8]!='F':
        pin+=1
    if cdr[9]!='S':
        pin+=1   
  
    if pin<max_pin:
        print('CDR_in pattern', cdr_in)
        # print('CDR_starts from', arr.index(cdr_in)+9)
        return (arr.index(cdr_in)+9)
 
    def cdr1_out(cdr_out):#WVRQAP
    
        pin=0
        max_pin=2            
        if cdr[1]!='V':
            pin+=1
        if cdr[0]!='W':
            pin+=1
        if cdr[2]!='R':
            pin+=1
        if cdr[3]!='Q':
            pin+=1    
        if cdr[4]!='A':
            pin+=1
        if cdr[5]!='P':
            pin+=1
            
        if pin<max_pin:
            # print('CDR_in pattern', cdr_out)
            # print('CDR_ends at', arr.index(cdr_out))
            return (arr.index(cdr_out))
 

K=10
arr=sequence
for i in range(len(arr)-k+1):
        slider=arr[i:k+i]
        print("CDR_1 is:", arr[cdr1_in(slider): cdr1_out(slider)])        
          

我是否正确地假设您正在分析免疫测序数据，并且 CDR 是指 B 或 T 细胞受体的互补决定区域？数据来自人类还是小鼠？如果是这种情况，您可能想看看现有的工具，而不是重新发明轮子。我用过mixcr https://github.com/milaboratory/mixcr。另一个流行的工具是IMGT/HighV-QUEST https://www.imgt.org/IMGTindex/IMGTHighV-QUEST.php但据我所知，它只能作为网络应用程序使用，不能用于大型数据集。如果它们不符合您的目的，您至少可能会得到有关如何继续的提示。