1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
|
/*****************************************************************************
*
* seqtest.c
* test program for sequence display, SeqPort and ToFasta
*
*****************************************************************************/
#include <seqport.h>
#include <tofasta.h>
void BuildList (SeqEntryPtr sep, Pointer data, Int4 index, Int2 indent);
#define NUMARG 2
Args myargs[NUMARG] = {
{"Filename for asn.1 input","stdin",NULL,NULL,FALSE,'i',ARG_FILE_IN,0.0,0,NULL},
{"Filename for output","stdout", NULL,NULL,TRUE,'o',ARG_FILE_OUT,0.0,0,NULL}};
Int2 Main(void)
{
AsnIoPtr aip;
SeqEntryPtr sep;
BioseqPtr PNTR seqlist;
Int4 seqnum, i, numseg, lens[10], j;
Int2 ctr;
SeqPortPtr spp;
Uint1 residue;
FILE* fp;
CharPtr title;
Char buffer[101];
MonitorPtr mon;
/* check command line arguments */
if ( ! GetArgs("SeqTest",NUMARG, myargs))
return 1;
mon = MonitorStrNew("SeqTest", 40);
SetProgMon(StdProgMon, (Pointer)mon);
/*
** Load SeqEntry object loader and sequence alphabets
*/
if (! SeqEntryLoad()) {
Message(MSG_ERROR, "SeqEntryLoad failed");
return 1;
}
/*
** Use the file "example.prt" as the ASN I/O stream. This file
** can be found in the ncbi/demo. It is in ASN.1 Print Value format.
*/
if ((aip = AsnIoOpen(myargs[0].strvalue, "r")) == NULL)
return 1;
/*
** Write the output to "seqtest.out".
*/
fp = FileOpen(myargs[1].strvalue, "w");
fprintf(fp, "Sequence summary:\n\n");
/*
** Read in the whole entry into the Sequence Entry Pointer, sep.
** Close the ASN stream, which in turn closes the input file.
*/
sep = SeqEntryAsnRead(aip, NULL);
aip = AsnIoClose(aip);
mon = MonitorFree(mon);
SetProgMon(NULL, NULL);
/*
** Determine how many Bioseqs are in this SeqEntry. Allocate
** enough memory to hold a list of pointer to all of these
** Bioseqs. Invoke an Explore function to "visit"each Bioseq.
** We are allowed to pass one pointer for use by the exploring
** function, in this case, "BuildList".
*/
seqnum = BioseqCount(sep);
seqlist = MemNew((size_t)(seqnum * sizeof(BioseqPtr)));
BioseqExplore(sep, (Pointer) seqlist, BuildList);
/*
** For each Bioseq in the SeqEntry write out it's title
** len, number of gaps, and number of segments. Write out
** the length of each segment, up to 10.
*/
for(i = 0; i < seqnum; i++) {
numseg = BioseqCountSegs(seqlist[i]);
title = BioseqGetTitle(seqlist[i]);
FilePuts((VoidPtr)title, fp);
FilePuts("\n", fp);
fprintf(fp, "len=%d gaps=%d segs=%d\n", BioseqGetLen(seqlist[i]),
BioseqGetGaps(seqlist[i]), numseg);
if ((numseg > 1) && (numseg <= 10)) {
BioseqGetSegLens (seqlist[i], lens);
for (j = 0; j < numseg; j++)
fprintf(fp, " len = %d\n", lens[j]);
}
FilePuts("\n", fp);
}
spp = SeqPortNew(seqlist[0], 0, -1, 0, Seq_code_iupacna);
if (spp == NULL)
Message(MSG_ERROR, "fail on SeqPortNew");
fprintf(fp, "SeqPort: plus strand with SeqPortGetResidue\n\n");
i = 0;
while ((residue = SeqPortGetResidue(spp)) != SEQPORT_EOF) {
if (! IS_residue(residue)) {
buffer[i] = '\0';
fprintf(fp, "%s\n", buffer);
i = 0;
switch (residue)
{
case SEQPORT_VIRT:
fprintf(fp, "[Gap]\n");
break;
case SEQPORT_EOS:
fprintf(fp, "[EOS]\n");
break;
default:
fprintf(fp, "[Invalid Residue]\n");
break;
}
}
else {
buffer[i] = residue;
i++;
if (i == 60) {
buffer[i] = '\0';
fprintf(fp, "%s\n", buffer);
i = 0;
}
}
}
if (i) {
buffer[i] = '\0';
fprintf(fp, "%s\n", buffer);
}
fprintf(fp, "[EOF]\n");
SeqPortFree(spp);
fprintf(fp, "\nSeqPort on minus with SeqPortRead\n\n");
spp = SeqPortNew(seqlist[0], 0, -1, Seq_strand_minus, Seq_code_iupacna);
if (spp == NULL)
Message(MSG_ERROR, "fail on SeqPortNew");
do {
ctr = SeqPortRead(spp, (Uint1Ptr)buffer, 60);
if (ctr > 0) {
buffer[ctr] = '\0';
fprintf(fp, "%s\n", buffer);
} else {
ctr *= -1;
switch (ctr)
{
case SEQPORT_VIRT:
fprintf(fp, "[Gap]\n");
break;
case SEQPORT_EOS:
fprintf(fp, "[EOS]\n");
break;
case SEQPORT_EOF:
fprintf(fp, "[EOF]\n");
break;
default:
fprintf(fp, "[Invalid Residue]\n");
break;
}
}
} while (ctr != SEQPORT_EOF);
SeqPortFree(spp);
/*
** Write out the nucleic acid sequences in this SeqEntry
*/
fprintf(fp, "\nNucleic Acids in FASTA format:\n\n");
SeqEntryToFasta(sep, fp, TRUE);
/*
** Write out the protein sequences in this SeqEntry.
*/
fprintf(fp, "\nProteins in FASTA format:\n\n");
SeqEntryToFasta(sep, fp, FALSE);
/*
** Close the output file and free up allocated space.
*/
fclose(fp);
MemFree(seqlist);
SeqEntryFree(sep);
return 0;
}
/*
** This SeqEntry exploration function copy the current pointer position inthe
** the Bioseq entry to a list of Bioseq pointers
*/
void BuildList(SeqEntryPtr sep, Pointer data, Int4 index, Int2 indent)
{
((BioseqPtr PNTR) data)[index] = (BioseqPtr)sep->data.ptrvalue;
return;
}
|