C语言读取fasta中核酸序列

原创 2016年05月30日 20:21:40

 背景介绍

在生物信息学中,FASTA格式(又称为Pearson格式),是一种基于文本用于表示核苷酸序列氨基酸序列的格式。在这种格式中碱基对或氨基酸用单个字母来编码,且允许在序列前添加序列名及注释。(来自百度百科)

其实fasta格式就是简单的文本格式,用C语言中可以用其提供的标准库函数对其进行操作,本文只简单介绍代码中所用到的函数;

文件的打开(fopen函数),函数原型:FILE * fopen(const char * path,const char * mode);打开文件实际是建立文件的联系,使文件指针指向该文件

文件的关闭(fclose函数),函数原型:int fclose( FILE *fp );关闭一个流,注意:使用fclose()函数就可以把缓冲区内最后剩余的数据输出到磁盘文件中,并释放文件指针和有关的缓冲区。

文件读(fread函数)函数原型:size_t fread(void* buff,size_t size,size_t count,FILE* stream);作用:从文件中读入数据到指定的地址中

文件重定位(rewind函数),函数原型:void rewind(FILE *stream); 作用:将文件内部的位置重新指向一个流数据流文件的开头

文件定位(fseek函数),函数原型:int fseek(FILE *stream, long offset, int fromwhere); 作用:定位流数据流文件上的文件内部位置指针

 程序功能介绍

本程序完成fasta文件的读取,并读取其中核酸序列以及一小段描述符,并将其地址存取于指针变量中,方便对其中每个字符进行操作,
下面贴上代码和注释。

#include <stdlib.h>
#include <stdio.h>
//#define NUM 10*1024
#define ERROR -1
char *buff;
char *start_point;
int size;
char *get_gi(char *p,int n){
	static char *m,*t;
	m=(char *)malloc(n);
	t=(char *)malloc(n);
	if(m==NULL||n==NULL)
		return NULL;
	t=m;
	while(n--){
	*m++=*p++;
	}
	*(m)='\0';
	return t;
}
char *get_seq(char *p,int n){
	static char *m,*t;
	m=(char *)malloc(n-1);
	t=(char *)malloc(n-1);
	if(m==NULL||n==NULL)
		return NULL;
	t=m;
	p=p+n;
	while((n<size)&&(*p)){
	*m++=*p++;
	n++;
	}
	*(m-1)='\0';
	return t;
}
int main(int argc, char* argv[])
{   
	FILE *file;
	size=0;
	file = fopen(argv[1],"rt");
	/*argv[1]为命令行传入的地址,"rt"表示以只读的方式打开一个文本文件*/
	if(file==NULL)
		{printf("ERROR:cannot retrieve this file.\n");
		return ERROR;}
	fseek(file,0L,SEEK_END);
	size=ftell(file);
	/*通过定位到文件末尾,读出文件大小size,或者也可通过下面注释掉的for循环读取文件大小size*/
	rewind(file);
	//char a;
	//a='0';
	//int ii;
	//ii=0;
	/*for(;a!=EOF;){
	a = getc(file);
	
	printf("a is %c%d\n",a,ii);
	ii++;
	}*/
	//rewind(file);
	//printf("ii number is %d\n",ii);
	printf("The file size is %d\n",size);
	buff = (char *)malloc(size-1);
	start_point = (char *)malloc(size-1);
	if(buff==NULL||start_point==NULL)
		return ERROR;
	fread(buff,size-1,1,file);
	/*将file指向的文本文件内容读入buff缓冲区中*/
	start_point=buff;
	/*start_point用于存储buff指向的首地址,用于free等*/
	printf("%s\n",buff);
	/*打印出文本文件内容,此处用于调试,printf是个很好的调试方法,此处可检查文本是否读出,以及是否正确等*/
	static int i;
	//unsigned short *aa;
	//printf("%c\n",*buff);
	static int pos;
	static int seq_pos;
	for(;*buff;buff++){
		//printf("%p\n",buff);
		i++;
		if((*buff=='|')&&(*(buff+1)==' ')){
			pos=i;
			//buff--;
			printf("The value of pos is %d\n",pos);
		}
		if((*buff=='A'||*buff=='T'||*buff=='C'||*buff=='G')&&(*(buff+1)=='A'||*(buff+1)=='T'||*(buff+1)=='C'||*(buff+1)=='G')\
			&&(*(buff+2)=='A'||*(buff+2)=='T'||*(buff+2)=='C'||*(buff+2)=='G')){
		seq_pos=i-1;
		printf("The value of seq_pos is %d\n",seq_pos);
		break;
		}
	}
	/*for循环中记录了标识符的结束位置和核酸序列的起始位置,这里的标识符是指的第一个空格前面的字符*/
	//printf("%d\n",i);
	char *mm=get_gi(start_point,pos);
	/*mm指向标识符的首地址,被调函数get_gi其实是用指针变量实现的字符串的拷贝,末了需要加上结束符'\0'*/
	char *seq=get_seq(start_point,seq_pos);
	/*seq指向核酸序列的首地址,被调函数get_seq也是指针变量实现的字符串的拷贝*/
	printf("DESCRIPTOR=\n%s\n",mm);
	printf("SEQ=\n%s\n",seq);
	buff=start_point;
	
	free(buff);
	//free(mm);
	//free(seq);
	fclose(file);
	//getchar();
	return 0;
}

编译执行

在Linux下编译通常通过makefile进行,它的一个好处就是自动化编译,文章此处不再赘述,网上关于makefile帖子非常多,我贴上这次用到的一个很简单的makefile ,源文件为seq.c,makefile注意和源文件在一个目录。

终端切换到当前目录下make,编译成功后会形成一个二进制可执行文件和中间目标文件seq.o

执行seq可执行文件,并传入要读取的文件

成功读出序列和标识符


结束语

C语言其实远没有想的简单,指针真的是C语言的灵魂所在,想要真正用好指针不容易。在python中读取这个序列很简单,就是调用模块的一个函数,但C语言总能让我们思考的更多,它不像python,java,标准库函数中可供我们调用的很有限,几乎每一步都要自己实现,很好地体现了面向过程的思想。

R/BioC序列处理之三:Biostrings模式匹配和序列比对

Biostrings最后一节,介绍模式匹配和序列比对的相关函数和操作。 下面我们使用拟南芥基因转录起始点上游1kb的序列进行分析。序列文件可以从TAIR网站(http://www.arabidopsi...
  • u014801157
  • u014801157
  • 2014年04月23日 16:21
  • 4052

Adam学习13之Fasta/Fastq/SAM/BAM文件格式数据读取

0.代码(读取方法): package org.bdgenomics.adamLocal.algorithms.test import org.apache.spark.SparkConf imp...
  • bob601450868
  • bob601450868
  • 2016年04月30日 22:33
  • 1912

R从文件中读取数据,输出文件

看了几天的书,终于到这一步了,说实话,用R来做统计,很少有人手动的去输入那些数字,肯定是从别的地方导入的,我们用来处理就可以了,所以到这里才算是真正的入门,前面都是做基础的练手。 我学习R从《R语言...
  • gaorongchao1990626
  • gaorongchao1990626
  • 2013年01月03日 16:35
  • 4819

读取DS18B20序列码的程序 C语言

  • 2015年11月14日 10:01
  • 30KB
  • 下载

核酸序列分析

转载一个很全面的核酸序列分析的文章 核酸序列分析 核酸序列分析 1、核酸序列检索 可通过NCBI使用Entrez系统进行检索,也可用EBI的SRS 服务器进行检索。在同时...
  • likelet
  • likelet
  • 2012年04月06日 15:03
  • 1755

从gff3文件获取fasta序列

chr1A NRGenome gene 1157233 1158291 . + . ID=TRIAE_CS42_U_TGACv1_641506_AA2096860.1.pa...
  • msw521sg
  • msw521sg
  • 2016年09月06日 11:13
  • 288

从gff3中获取fasta序列

#!/usr/bin/env python # -*- coding: utf-8 -*-from Bio import SeqIO record_dict = SeqIO.index(open('g...
  • msw521sg
  • msw521sg
  • 2016年09月05日 23:11
  • 263

FASTA序列格式说明

fasta序列格式是blast组织数据的基本格式,无论是数据库还是查询序列,大多数情况都使用fasta序列格式,所以首先对fasta格式在做详细说明。 下面是一个来源于NCBI的fasta格式序列: ...
  • SHMILYRINGPULL
  • SHMILYRINGPULL
  • 2011年12月26日 10:04
  • 3012

matlab写数据到txt文件,C语言读取文件内容到数组@项目简介:基于PSS序列(频域)估计整数倍频偏

一、编程思想: 1.以文本方式打开文件。 2.循环用fscanf格式化输入数据到数组。 3.判断fscanf的返回值,如果显示到达文件结尾,退出输入。 4.关闭文件。 5.使用数据。 二、代码实现:...
  • GSH_Hello_World
  • GSH_Hello_World
  • 2016年07月26日 11:20
  • 995

system verilog与c语言接口读取yuv图像序列

system verilog代码 `define PIC_W 176 //图像宽 `define PIC_H 144 //图像高 `define FRAME 300 //总帧数 `define ...
  • qin_lin_sb
  • qin_lin_sb
  • 2012年10月14日 23:40
  • 1656
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:C语言读取fasta中核酸序列
举报原因:
原因补充:

(最多只允许输入30个字)