使用位置文件和目标文件提取特定位置的序列

Question

我有一个 DNA 序列文件 1（250M 字符/字节），看起来像这样（FASTA 格式）：

$sequence-file1
TCCTCCAAATGATGTCAGTGTCCTCCATATGATGTCAATGTCCTCCATAT
GATGTCAATATCCTCCGTATGATGTCAATATCCTCCGTATGATGTCAATA
TCCTCCATATGATGTCAGTGTCCTCTGTATGACATCAATATCCTCCATAC
GATGCCCCTGTCCTTCATATGATGTCAGTGTCCTTTTGTGAGCACCAGTG
TCCTTTGTATGACATCAGTAGTCTCCCATGAATGTCACTGTCTTCCCATA

以及这种格式的序列位置文件2，位置不连续：

$positions-file2
1
2
7
39
51

我需要从sequence-file1中提取positions-file2中指定位置的字符并打印出“位置字符”如下awk程序：

$prog.file.awk    
    {
        for (i=1;i<=length;i++) 
            if((i+(NR-1)*length)==x) 
                print x"\t"substr($0,i,1);exit 
    }

...当我通过xargs 将“x”的位置传递给它时，仅对前 50 行执行此操作： xargs -I{i} awk -v x={i} -f prog.file.awk sequence-file1 < positions-file2 输出：

1   T
2   C
7   A
39  T

positions-file2 中任何大于 50 的数字都将被忽略。给定上述输入文件，我想要的输出是：

1   T
2   C
7   A
39  T
51  G

我也在寻找一种经济的解决方案，因为对于 250M 的字符文件，我有大约 200M 的位置要匹配。

Answer 1

在您更正数据后，以下操作将起作用：

awk 'FNR==NR{p[$1]++;next} {for(x in p)print x,substr($0,x,1)}' pf2 sf1

目前每行只有 50 个字符，因此您无法打印第 51 个字符。它也不会搜索行中的每个字符，它只是提取您指定的字符，因此会快得多。

说明

FNR==NR 表示后面花括号中的所有内容仅适用于文件pf2 的处理。在那里，我将位置保存在数组p[] 中，因此在读取位置文件后 p[1]=1, p[2]=1, p[7]=1, p[39]=1 and p [51]=1。

第二组花括号中的代码仅适用于第二个文件sf1。它循环遍历我们保存在p[] 中的所有位置，并通过使用substr() 提取它们来打印当前记录中选定的字符。

Answer 2

我知道标签上写着 awk，但考虑到数据集的预期大小，awk 感觉像是错误的工具。我的 C 最终比预期的要长一点，但部分是因为我添加了代码来验证行终止和行长。

[dennis@localhost dna]$ gcc -Wall reindex.c 
[dennis@localhost dna]$ ./a.out sequence.dat position.dat
1   T
2   C
7   A
39  T
51  G

在我将您的 sequence-file1 示例文本复制到 sequence.dat 和 position-file2 到 position.dat。

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <errno.h>


void usage(int argc,char **argv);
int analyze(
  FILE *fp 
  ,long *pLineTextLen   /**< OUT: Length of alpha text per line      */
  ,long *pLineBinLen    /**< OUT: Total length of line including lf  */
  );

int reindex(
  FILE *seqFp           /**< IN: file with sequence to reindex       */
  ,FILE *posFp          /**< IN: file with indexes to extract        */
  ,long lineTextLen     /**< IN: text to index per line              */
  ,long lineBinLen      /**< IN: characters including termintion     */
  );


int main( int argc, char **argv)
{
  int errval;
  FILE * seqFp=NULL;
  FILE * posFp=NULL;
  long   lineTextLen;
  long   lineBinLen;
  char  *sequenceName=NULL;
  int    argIdx;

  argIdx=1;

  if(argIdx >= argc)
  {
    usage(argc,argv);
    errval=-__LINE__;
    goto exiterror;
  }
  seqFp = fopen(argv[argIdx],"r");
  if(seqFp == NULL)
  {
    errval=errno;
    fprintf(stderr,"Unable to open %s\n",argv[argIdx]);
    goto exiterror;
  }
  sequenceName = argv[argIdx];
  argIdx++;
  if(argIdx >= argc)
  {
    usage(argc,argv);
    errval=-__LINE__;
    goto exiterror;
  }
  posFp = fopen(argv[argIdx],"r");
  if(posFp == NULL)
  {
    errval=errno;
    fprintf(stderr,"Unable to open %s\n",argv[argIdx]);
    goto exiterror;
  } 
  errval = analyze(seqFp,&lineTextLen,&lineBinLen);
  if(errval)
  {
    fprintf(stderr,"Unable to estimate line length of %s\n"
            ,sequenceName);
    errval=-__LINE__;
    goto exiterror;
  }
  errval = reindex(seqFp,posFp,lineTextLen,lineBinLen);
  if(errval)
  {
    fprintf(stderr,"Unable to reindex (errval=%i)\n"
            ,errval);
    goto exiterror;
  }


exiterror:
  if(seqFp != NULL)
  {
    fclose(seqFp);
    seqFp=NULL;
  }
  if(posFp != NULL)
  {
    fclose(posFp);
    posFp=NULL;
  }
  return(errval);

}


void usage(int argc,char **argv)
{
  (void)argc;  /* yes I'm ignoring it atm */

  fprintf(stderr,"%s {seqeuence-file} {position-file}\n"
          ,argv[0]);
  return;
}

/*********************************************************************/
/** Analyze file to determine line lenth
 * 
 * Analyze first few lines of file for identical length text
 * lines consisting only of alpha text.
 * 
 * return non-zero if lines not consistent or other error.
 *********************************************************************/
int analyze(
  FILE *fp 
  ,long *pLineTextLen   /**< OUT: Length of alpha text per line      */
  ,long *pLineBinLen    /**< OUT: Total length of line including lf  */
  )
{
  int input;
  int lineTextLen=0;
  int lineBinLen=0;
  int confirmCount=0;
  int count=0;
  enum
  {
    TEXT_READ=0,
    TERM_READ=1
  }
  state= TEXT_READ;

  do
  {
    input=fgetc(fp);
    if(input != EOF)
    {
      if(isalpha(input))
      {
        if( state == TERM_READ)
        {
          state = TEXT_READ;
          if(lineBinLen != 0 )
          {
            if( count != lineBinLen )
            {
              /* mismatch */
              goto exiterror;
            }
            confirmCount++;
          }else
          {
            lineBinLen=count;
          }
          count=0;  /* start new line */
        }
        count++;
      }
      else if( ( input == '\r' )
               || (input == '\n')
               || isblank(input) )
      {
        if(state == TEXT_READ)
        {
          state = TERM_READ;
          if(lineTextLen!=0)
          {
            if(lineTextLen  != count )
            {
              /* mismatch */
              goto exiterror;
            }
            confirmCount++;
          }
          else
          {
            lineTextLen=count;           
          }
        }
        count++;
      }
    }
  }
  while(input!=EOF 
        && confirmCount<4); /* 2 text and 2 bin */
exiterror:  
  rewind(fp);
  if( pLineTextLen )
  {
    *pLineTextLen = lineTextLen;
  }
  if( pLineBinLen )
  {
    *pLineBinLen = lineBinLen;
  }

  return(confirmCount<4);  /* non-zero if not confirmed */
}

/**********************************************************************/
/** reindex sequence file to std out.
 * 
 * Print char at specified character indexes in sequence file.
 * Character indexes are one-based index of characters in 
 * seq file not including line terminations.  Line length and 
 * termination are assumed to be consistent and specified by 
 * passed parameters.
 *
 * Indexes are read as text strings one per line from pos file.
 *
 * /return non-zero on error. 
 *********************************************************************/

int reindex(
  FILE *seqFp           /**< IN: file with sequence to reindex       */
  ,FILE *posFp          /**< IN: file with indexes to extract        */
  ,long lineTextLen     /**< IN: text to index per line              */
  ,long lineBinLen      /**< IN: characters including termintion     */
  )
{
  int  errval=0;
  char buffer[80];
  char *pInput=NULL;
  long  index;
  long  lines;
  long  seekPos;
  int   sequence;

  do
  {
    pInput=fgets(buffer,sizeof(buffer),posFp);
    if( (pInput != NULL)
        && ( !isalnum(pInput[0]) ))  /* empty line */
    {
      pInput=NULL;
    }

    if(pInput != NULL)
    {
      index=strtol(pInput,NULL,0);
      if(index==0)
      {
        errval=-__LINE__;
        goto exiterror;
      }
      index--;  /* switch to zero based index */
      /* integer truncated division expected below */
      lines=index/lineTextLen;
      seekPos= ( ( lines * lineBinLen ) 
                 + ( index - lines * lineTextLen ) );

      fseek(seqFp,seekPos,SEEK_SET);
      sequence=fgetc(seqFp);
      if(sequence == EOF)
      {
        errval=-__LINE__;
        goto exiterror;
      }
      fprintf(stdout,"%li\t%c\n"
              ,index+1 /* convert back to one based */
              ,sequence);
    }

  }while(pInput!=NULL);

exiterror:
  return(errval);

}