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,标准库函数中可供我们调用的很有限,几乎每一步都要自己实现,很好地体现了面向过程的思想。

版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

核酸序列分析

转载一个很全面的核酸序列分析的文章 核酸序列分析 核酸序列分析 1、核酸序列检索 可通过NCBI使用Entrez系统进行检索,也可用EBI的SRS 服务器进行检索。在同时...

核酸序列分析

核酸序列分析 核酸序列分析 1、核酸序列检索 可通过NCBI使用Entrez系统进行检索,也可用EBI的SRS 服务器进行检索。在同时检索多条序列时,可通过罗逻辑关系式按照GenBank接受号进...

FASTA序列格式说明

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

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

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

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

system verilog代码 `define PIC_W 176 //图像宽 `define PIC_H 144 //图像高 `define FRAME 300 //总帧数 `define ...

从gff3文件获得fasta序列

get fasta sequence from gff3 chr1A NRGenome exon 1157233 1158291 100 + . ID=TRIAE_CS42_U_TGACv1_641...

生物信息脚本练习(1) 找出fasta文件中大于500的序列

最近做了一些生物信息的脚本练习。 这是第一个例子。 找出一个fasta文件中大于500的序列,并重定向到另一个新的文件中。 这个文件每条序列是如下的样子。 c100027.graph_c0|...
  • hxoxh
  • hxoxh
  • 2017-07-30 15:38
  • 138

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

0.代码(读取方法): package org.bdgenomics.adamLocal.algorithms.test import org.apache.spark.SparkConf imp...

C语言序列点总结

  • 2014-06-19 10:25
  • 1.25MB
  • 下载
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:深度学习:神经网络中的前向传播和反向传播算法推导
举报原因:
原因补充:

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