python - 在 NGS 数据中寻找 CDR

Question

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

序列:-'FYSHSAVTLDESGGGLQTPGGGLSLVCKASGFTFSSYGMMWVRQAPGKGLEYVAGIRNDA GDKRYGSAVQGRATISRDNGQSTVRLQLNNLRAEDTGTYFCAKESGCYWDSTHCIDAWGH GTEVIVSTGG'。

cdr1 从:- 'VCKASGFTFS' 开始，最多可替换三个，但 C 必须排在第二位。cdr1 结束于:-'WVRQAP'，最多有两个替代品，但 R 必须在第 3 位。

提取的 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)])        
          

Answer 1

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