cannin/biostars_qa huggingface dataset:
Question: IGB has previously worked on my ubuntu 15.04 desktop, but within the last week it will load a genome from the default options, but die upon trying to load another track. Sometimes it will just crash during the time between loading the genome and doing anything else. I've tried with bed, gff3 and wig files. It appears that it has loaded the file from the log, but then crashes with the following error:

```
#
# A fatal error has been detected by the Java Runtime Environment:
#
# SIGSEGV (0xb) at pc=0x00007f3de2676714, pid=12777, tid=139903480776448
#
# JRE version: OpenJDK Runtime Environment (8.0_45-b14) (build 1.8.0_45-internal-b14)
# Java VM: OpenJDK 64-Bit Server VM (25.45-b02 mixed mode linux-amd64 compressed oops)
# Problematic frame:
# C [libawt_xawt.so+0x45714]
#
# Failed to write core dump. Core dumps have been disabled. To enable core dumping, try "ulimit -c unlimited" before starting Java again
#
# An error report file with more information is saved as:
# /home/sbi6dap/Local/IGB/hs_err_pid12777.log
#
# If you would like to submit a bug report, please visit:
# http://bugreport.java.com/bugreport/crash.jsp
#
```

My java is 1.8 as is required (see above).

Has anyone else had this difficulty?

Edit: I've tried increasing the allocated ram with the `-Xmx4g` option in the `.vmoptions` file, but also turned it back to `-Xmx1280m` in consideration of 32bit limits (although I am running 64). My machine has 32gb available.

Answer: Hi Daniel,

I'm very sorry to hear you are dealing with this bug. This is a problem specific to some Ubuntu versions which we have now resolved on our development branch. It is possible to work around this bug by dragging and dropping files to open them instead of using the normal "open file" action; however, we will be releasing an updated version of IGB very soon (early next week time frame) with a patch.
Question: I tried to fit a linear model on a set of metabolomics data for an experiment: drug treatment for a group of animals to receive 7 different types of concentration. Each group of animals there are 10 replicates. I would like to plot the changes of a particular metabolites across different concentration of drug given. It is like a factorial design here.

Everything works well until it comes to the end when I wanted to get the logFC which represents the gradient of the overall changes.

grp=factor(target$Group, levels=c("0","0.1","10","30","50","80","100"))
design <- model.matrix(~grp)
fit <- lmFit(data, design)
topTable(fit, n = Inf, adjust="BH")

The output gave me a list of coefficients which corresponds to each level of the design , a column of AveExpr, F, P.value and adj.Pval.
I wonder if anyone can kindly suggest me a way to get the logFC? either calculate manually or is there any changes on the code that I have missed out? I saw many posts mentioned that limma is treating each time point independently in a time series experiment.
I am not interested in contrast and treating each group of treatment separately.
I wonder if is there a way for me to calculate the logFC in this case.

The set up of the experiment is somehow similar to the post below, but in a simpler way
https://www.biostars.org/p/157074/

Answer: limma has computed logFCs for you. It has computed logFC for dose 0.1 vs 0, for dose 10 vs 0, for 30 vs 0 and so on. That is how time course data is usually analysed. The logFC of dose 100 vs dose 0 is a sort of overall logFC because it measures the effect of the highest dose vs the lowest.

Your email sounds that like you instead wanted to fit a linear trend to the log expression values as a function of dose (or of log dose perhaps). That is easy to do in limma but it is not good practice to assume that log expression is a straight line function of dose (or of log dose perhaps). I advise against it.
Question: I am trying to align the following sequence using a python code. I am able to obtain the best score results but not the alignments with lower scores. Is there an efficient way to do this?

# Import pairwise2 module
from Bio import pairwise2

# Import format_alignment method
from Bio.pairwise2 import format_alignment

# Define two sequences to be aligned
X = "ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGGATCGGAATATCACTATTCTGCTCGAACTTCCATAG"
Y = "ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGAATATCACTATTCTGCTCGAACTTCCATAGAGATCG"

# Get a list of the global alignments between the two sequences ACGGGT and ACG satisfying the given scoring
# A match score is the score of identical chars, else mismatch score.
# Same open and extend gap penalties for both sequences.
alignments = pairwise2.align.globalms(X, Y, 1, -1, -1, -0.1, force_generic=True, score_only = True)

# Use format_alignment method to format the alignments in the list
for a in alignments:
print(format_alignment(*a))

Here is the result I get with the above code.

ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGGATCGGAATATCACTATTCTGCTCGAACTTCCATA------G
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||| |
ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGT------GAATATCACTATTCTGCTCGAACTTCCATAGAGATCG
Score=141

What I want is alignments with lower scores such as 130, 100, etc. I know there could be many such examples but I need the alignments with the lower scores. I tried looking at the source code for Biopython but it seems to only give the best score.

Answer: In case it wasn't clear what I meant above, here is how making artificial gaps would look. Real alignment first:

ACAGCCAGCGGGGGAAATAGA
|||||||||||||||||||||
ACAGCCAGCGGGGGAAATAGA

Now the alignment with one artificial move in up and left directions during the backtracking:

ACAG-CCAGCGGGGGAAATAGA
||| |||||||||||||||||
ACA-GCCAGCGGGGGAAATAGA

This alignment will have a worse score, but it should also be clear that it is meaningless. Depending on the length of the two sequences, one can create literally millions of these substandard and meaningless alignments, especially if the gap is allowed to be longer.

This is different from legitimately different alignments of this kind, which all have the same score:

AGGGCT
| ||||
A-GGCA

AGGGCT
|| |||
AG-GCA

AGGGCT
||| ||
AGG-CA
Question: Hi,

I have bacterial RNA-seq data from paired-end reads (2x75). I'm interested in sense/antisense differential expression.

I would like to know if makes sense to merge R1 and R2 into one read (if overlaps) or maybe work only with 75 bp reads (R2 ones).

Thanks, Bernardo

Answer: Why bother? Presuming that you used a stranded protocol when making the libraries, you can just align things as pairs and tell whichever program you use to do counting to observe strandedness of pairs. Both htseq-count and featureCounts can do that.

To get at the antisense transcription, one could just run the counters again and just tell them to use the other read as the strand reference (e.g., running htseq-count -s yes ... and then htseq-count -s reverse ...).
Question: Kind of a basic awk question but just trying to take the following intervals:

chr1 10001 10150
chr1 10008 10090
chr1 10025 10207
chr1 10058 10172

And generate a file output with a single command that in the first column has the midpoint of each and in the second column next to it has the entire interval length, such as this:

10075 149
10049 82
10116 182
10115 114

Thanks.

Answer: $ awk '{print int(($2+$3)/2), $3-$2}' test.txt
10075 149
10049 82
10116 182
10115 114
Question: Dear all,

Do you have any idea how to easy extract contigs from fasta file wich contains specific sequence?

For example:

My sequence:

ACCGTACCC

My FASTA:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
>c1044
ATGAATAAAATAATTTTGTATCATAAATCGAGCTGTTAATTATT
>c1044
TTCATATTTGTAGCTAAGCAGAGGCGAAGCGTTCTTGTATCG
```

My output:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
```

Thank you so much for any ideas and help.

Answer: You can use `sed` + `grep` (as suggested by NicoBxl and Devon) or BioPerl/BioPython (as suggested by Devon) or [Heng Li's bioawk][1]:

bioawk -c fastx '$seq ~ /ACCGTACCC/ { print ">"$name"\n"$seq; }' #might need a bit of tweaking

[1]: https://github.com/lh3/bioawk
Question: I am attempting to use TopHat2 (and, consequently, Bowtie2) on a computer cluster. Since it's on someone else's cluster, I did not do the installation of Bowtie, nor do I have control of the bowtie2-index.loc files.

I have a directory, .../Assemblies/hg19/, containing these files (the index files/reference files for Bowtie2):

```
hg19.1.bt2
hg19.3.bt2
hg19.2.bt2
hg19.4.bt2
hg19.rev.1.bt2
hg19.rev.2.bt2
hg19.zip
make_hg19.sh
```

I point to these index files using `tophat2 -G [indexFileDirectory]`. However, I cannot get tophat2 to register the files appropriately.

For instance, if I write `[indexFileDirectory]` as:

.../Assemblies/hg19/hg19

(The `...` represents a directory path), then I get this error:

Error: cannot find transcript file .../Assemblies/hg19/hg19

Here's a table of corresponding index file location and errors

```
[indexfileDirectory] Error
…/Assemblies/hg19/ Error: is a directory
…/Assemblies/hg19/hg19 Error: cannot find transcript file .../Assemblies/hg19/hg19
…/Assemblies/hg19/* Error: Could not find Bowtie 2 index files (.../Assemblies/hg19/hg19.2.bt2.*.bt2)
…/Assemblies/hg19/*hg19* Error: Could not find Bowtie 2 index files (.../Assemblies/hg19/hg19.2.bt2.*.bt2)
…/Assemblies/hg19/hg19* Error: Could not find Bowtie 2 index files (.../Assemblies/hg19/hg19.2.bt2.*.bt2)
```

My review of internet forums sounds like the second option (`.../hg19/hg19`) should be the one that works. I'm going nuts! Did anyone else have this problem?

Answer: If you see tophat help:

```
Usage:
tophat [options]
```

You are trying to put the bowtie index with `-G` parameter, and `-G` is for specifying an annotation file.

```
-G/--GTF (GTF/GFF with known transcripts)
```
Question: Hey. Is there quick way (biopython or something) to retrive secondary structure from dssp file?

eg 'HHHHHHHHHHEEEEEEEEEEEECCCCCCCCHHHHHHEEEEEEEEE'

I need to get SS by pdb id and build FASTA like records

>3gja
HHHHHHHHHHEEEEEEEEEEEECCCCCCCCHHHHHHEEEEEEEEE
>3akl
CCCCCCCCCCCCEEEEEEEEEEEHHHHHHHHHHEEEEEEEEEEEECCCCCCCCHHHHHHEEEEEEEEE
...

Answer: Here is Biopython's [**DSSP code**][1], with typical use shown near the top. It is not difficult to extract this information on your own by reading contents of column 17 of the DSSP output, following this line:

# RESIDUE AA STRUCTURE BP1 BP2 ACC N-H-->O O-->H-N N-H-->O O-->H-N TCO KAPPA ALPHA PHI PSI X-CA Y-CA Z-CA

You may want to convert the 8-letter designation utilized by DSSP into 3 letters - space (or empty character `''`) is considered to be `C`. Typical conversion schemes are:

SS-Scheme 1: H,G,I->H ; E,B->E ; T,S->C
SS-Scheme 2: H,G->H ; E,B->E ; I,T,S->C I think this is most common
SS-Scheme 3: H,G->H ; E->E ; I,B,T,S->C
SS-Scheme 4: H->H ; E,B->E ; G,I,T,S->C
SS-Scheme 5: H->H ; E->E ; G,I,B,T,S->C

[1]: https://github.com/biopython/biopython/blob/master/Bio/PDB/DSSP.py
Question: I want to publish a script for general public use under GPL licensing. My question is how should I go about it?

I googled about the same and my understanding is:

1. I should paste the following text in comment section of my script:

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see .

2. I should provide a license.txt file containing GPL license terms and conditions along with my script copied from following link:
https://www.gnu.org/licenses/gpl-3.0.txt

Would these 2 points sufficient? Or do I need to contact somebody from GPL?

Answer: Actually, you need to have an additional first paragraph (see [How to use GNU licenses for your own software][1]):

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see .

That's it. There is no need to contact anyone or to ask for permission.

For single-file-scripts, I personally do not attach the license itself but modify the last paragraph to something like:

You can receive a copy of the GNU General Public License at .

[1]: https://www.gnu.org/licenses/gpl-howto.html
Question: Is there a way to only pass certain fields from a `record` type to a step? As opposed to the whole record object?

This is my best guess, but I was wondering if there was a cleaner way:

inputs:
samples:
type:
type: array
items:
type: record
fields:
fastq1: File
fastq2: File
SM: string

steps:
step_1:
in:
samples: samples
out: [bams]

step_2:
in:
bams: step_1/bams
SM:
valueFrom: ${ return inputs.samples.map(function(x) { return x.SM }) }
out: [changed_bams]

Answer: You almost have it:

steps:
step_2:
in:
SM:
source: samples
valueFrom: ${ return self.map(function(x) { return x.SM })
Question: Although I already found some ways to annotate genomic positions with rsIDs using e.g. UCSC table browser, I'm not happy with that since I want a one-in-all linux script taking also strand issues (flipped alleles A-T vs- T-A or switched reference alleles) into account.

What I have:

chr position ref alt
10 169560 G T
10 171117 G A
10 171126 G A
10 172995 A C
10 178499 C T

What I want:

chr position ref alt rsID
10 169560 G T rsXXX
10 171117 G A rsXXX, rsXXX
10 171126 G A rsXXX
10 172995 A C rsXXX
10 178499 C T rsXXX

Thanks

Answer: I will write down my solution as an answer for documentation purposes. I started as Pirerre recommended, but then I used bcftools instead of GATK.

First, I created a header .txt file for the custom vcf file

##fileformat=VCFv4.0
##fileDate=09052019
##source=allchr_allvsall_sex_adjusted
##reference==GRCh37.p13
##phasing=partial
##INFO=
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO

Then I used `awk` to generate the data for vcf according the specifications (8 columns). Setting ID="." == missing, Quality to 100 and PASS for the filter for all positions. Of note my_chr_pos_alt_ref.out.gz data consists only of autosomal SNVs!

zcat my_chr_pos_alt_ref.out.gz | awk '{print $1, ".", $2, $3, $4, 100, "PASS", "AA="$3}' OFS='\t' > tmp.vcf

add the header

cat header.txt tmp.vcf > mydata.vcf
rm tmp*

zipped and indexed

bgzip mydata.vcf
tabix -p vcf mydata.vcf.gz

Finally annotated rsIDs using:

bcftools annotate \
-a 00-common_all.vcf.gz \
-c ID mydata.vcf.gz \
--output-type z \
-o mydata_dbSNP151.vcf.gz

dbSNP files from ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/
Question: I am a computer engineering student currently doing a research about the Trinity platform. I've chosen arabidopsis thaliana as my model organism. I downloaded small sets of data from SRA to test, but I need a large amout of raw data for my purpose. From where can I download all the raw transcriptome data for the species as a whole?

Answer: You could try an Entrez query like http://www.ncbi.nlm.nih.gov/sra?term=((Arabidopsis%20thaliana%5BOrganism%5D)%20AND%20%22biomol%20rna%22%5BProperties%5D)%20AND%20%22study%20type%20transcriptome%20analysis%22%5BProperties%5D . This produces a list of almost 5000 Arabidopsis transcriptomes accessions. On the top of your Entrez result page you can send the results to the "Run Selector" and then download the list of all Arabidopsis transcriptomes accessions and feed that into the SRA toolkit. It looks like about 6 Tb of data in all.
Question: Hey
I already did automatic genome annotation using RAST software. I am looking for simple visualisation tool or anything which can provide me information about ORF total number of given dataset or would let me extract decription of all ORFs together.

It's pretty simple to do, but programs I am using are not counting orfs and they won't let me copy all ORFs description where I could just number it on word or any other tool. I am using genebank format.

I am not looking for tools like ORF FINDER, orfs are already predicted and translated in my dataset.

Answer: If you have a FASTA file for your proteins, this is a simple way to count proteins:

grep ">" your_file_name.fas | wc -l

To extract annotations and save them in `annotations.txt`:

grep ">" your_file_name.fas > annotations.txt

If you want to annotate your genome by metabolic or onthology criteria, [**KEGG**][1] will do the trick.It can give you an output that looks like this:

![enter image description here][2]

[1]: https://www.genome.jp/kegg/
[2]: https://ncgas.org/Blog_Posts/Blog%20Images/Complete.png
Question: I'm just a beginner in this area so I apologize first for wasting your time with this stupid question:

I recently downloaded some data from Uniprot and want to do some domain coverage analysis, but just realized that there are some overlapping ranges that needs to be merged.

For example, my data looks like this:

no.	pfamseq_acc	seq_version	crc64	md5	pfamA_acc	seq_start	seq_end
1857719	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	11	244
1857720	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	482	737
1857722	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	83	299
1857723	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	286	500
1857724	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	10	222
1857725	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	223	423

And what I want is:

no.	pfamseq_acc	seq_version	crc64	md5	pfamA_acc	seq_start	seq_end
1857719	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	11	244
1857720	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	482	737
1857722	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	83	500
1857724	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	10	423

I have tried to write it in R which looks like:

library(dplyr)
setwd("C:/Users/hy/Documents/Desktop/test2")
d<- read.csv("test2.tsv", header = T, sep = "\t", dec = ".")

d %>%
group_by(md5) %>%
arrange(md5,seq_start) %>%
mutate(indx = c(0, cumsum(as.numeric(lead(seq_start)) >
cummax(as.numeric(seq_end)))[-n()])) %>%
group_by(md5,pfamseq_acc,pfamA_acc)
summairse(start = min(seq_start), end = max(seq_end))

The "summarise" could give me the correct result, but I failed to write it in a new file.

So I would appreciate it if anyone could give some ideas about how to achieve that (preferably by R or Python), any other methods is also welcomed.

Thanks.

Answer: Using Unix tools and BEDtools, reading a file `in.tsv`:

## Write coordinates as BED with the remaining information as one string in column 4
awk 'OFS="\t" {if ($1 == "no.") next} {print $2, $7, $8, $1"--"$6"--"$3"--"$4"--"$5}' in.tsv | \

## sort
sort -k1,1 -k2,2n -k3,3n | \

## save intermediate file for later intersection
tee tmp_sorted.tsv | \

## merge overlapping intervals
bedtools merge -i - | \

## intersect the merged interval with the intermediate file to restore additional information
bedtools intersect -a - -b tmp_sorted.tsv -wa -wb | \

## cut out what is relevant
cut -f1-3,7 | \

## get unique list
sort -k1,1 -k2,2n -u | \

## write back to original format
awk 'OFS="\t" {gsub("--","\t");print $4, $5, $6, $7, $8, $1, $2, $3}' | \

## add header
cat <(echo 'no. pfamseq_acc seq_version crc64 md5 pfamA_acc seq_start seq_end' | tr " " "\t") /dev/stdin > out.tsv && rm tmp_sorted.tsv
Question: Hi all,

I'm trying to convert minimal dosages into beds using gcta. Here's my code:

gcta64 --dosage-mach chr22.imputed.dose chr22.imputed.info --make-bed --out chr22.fdgsdf

My error is:

Reading dosage data failed. Please check the format of the map file.

Here's what my info file looks like:

```
SNP REF ALT Major Minor MAF AvgCall Rsq Genotyped LooRsq EmpR EmpRsq Dose1 Dose2
22:16050115 G A G A 0.00639 0.99361 0.00000 - - - - - -
22:16050213 C T C T 0.00759 0.99241 0.00000 - - - - - -
22:16050568 C A C A 0.00040 0.99960 0.00005 - - - - - -
22:16050607 G A G A 0.00100 0.99900 0.00002 - - - - - -
22:16050627 G T G T 0.00040 0.99960 0.00005 - - - - - -
22:16050654 A A 0.00180 0.99820 0.00000 - - - - - -
22:16050654 A A 0.01737 0.98263 0.00000 - - - - - -
```

Any ideas?

Best

Answer: Answer:

```
for chr in `seq 1 22`;
do
awk -v OFS='\t' '!($4="") , !($5="")' chr$chr.imputed.info > chr$chr.imputed.info.reworked
sed -i 's/\t\t//' chr$chr.imputed.info.reworked
awk '{t=$5; $5=$4; $4=t; print}' chr$chr.imputed.info.reworked > t && mv t chr$chr.imputed.info.reworked
awk '{$6=.999 FS $6;}1' OFS='\t' chr$chr.imputed.info.reworked > t && mv t chr$chr.imputed.info.reworked
awk '$1=$1"_"$2"_"$3"_"++i' chr$chr.imputed.info.reworked > t && mv t chr$chr.imputed.info.reworked
sed -i '1 s/^.*$/SNP\tAl1\tAl2\tFreq1\tMAF\tAvgCall\tRsq\tGenotyped\tLooRsq\tEmpR\tEmpRsq\tDose1\tDose2/g' chr$chr.imputed.info.reworked

gzip chr$chr.imputed.info.reworked
/labseq/tools/gcta_1.24.2/gcta64 --dosage-mach-gz chr$chr.imputed.dose.gz chr$chr.imputed.info.reworked.gz --make-bed --out chr$chr.imputed
done
```
Question: Hello,

I have a question regarding the RNA-star alignment program. I'm trying to map thousands of reads (~7000 sequences) against a fasta file containing a handful of sequences (~20 sequences). In order to do that in RNA-star, I had to generate a genome index based on the small fasta file that contains ~20 sequences:

$STAR --runMode genomeGenerate --genomeDir $direc --genomeFastaFiles $filename --runThreadN 4

Next I run the following command to do the mapping:

$STAR --genomeLoad NoSharedMemory --outFileNamePrefix $direc/STAR. --genomeDir $direc --readFilesIn file.fasta

After running the above I get the following segmentation fault error:

```
Jul 22 11:44:49 ..... Started STAR run
Jul 22 11:44:50 ..... Started mapping
Segmentation fault (core dumped)
```

I think this has something to do with the size of the genome index I'm generating. This is because I was able to run the above commands successfully when I had a larger number of sequences used to generate the genome index. Has anyone run into this issue before? If yes, any ideas/suggestions on how to resolve it?

Answer: I found what was the problem. When generating a genome index, this parameter must be changed to accommodate smaller genomes:

genomeSAindexNbases 14

> int: length (bases) of the SA pre-indexing string. Typically between 10 and 15. Longer strings will use much more memory, but allow faster searches.
Question: Hi.

I have a fasta file with multiple headers and want to get reverse complement sequences.
I usually use FASTX-TOOLKIT, but I want to learn how to do with linux commands.

I tried it using 'awk' and 'if' something like below, but I got a tedious result.

awk '{if($1 ~ />/) print $0; else gsub("[ATCG]","[TAGC]");print $0}' TEST.fa

As you can see, the file has the same string in the header.

How can I get reverse complement sequences without altering the headers?

> ATCG_1
agtcgacATCGATtataggta
> ATCG_2
cgactgcagtcgacATCGACT

Thank you!

Answer: using '**tr**', and '**rev**', 2 lines per sequence

cat input.fa | while read L; do echo $L; read L; echo "$L" | rev | tr "ATGC" "TACG" ; done