Hot answers tagged

15

Thanks to Manu Tamminen for this solution: echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev


8

All you need is cat. You won't find any better tool for a simple job like this. Just run: cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz cat NA24694_GCCAAT_L001_R4*fastq.gz > ...


8

awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa An explanation of the non-obvious bits: BEGIN{FS="::"} split columns using :: as the designator $1~">" If there's a > in the line substr($1,2) Trim the > that was at the beginning of the line.


8

For a quick estimate you’re making it more complicated than necessary. The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size. samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list. Putting this together, ...


7

I think it is easy to do with awk, e.g. see here. For your example, with your data in file: % awk '$2 != $3 {print $0}' file # outputs cont100 1128 1125 cont1006 103 19 cont104 895 890 cont1078 43 42 cont1081 727 2 % awk '$2 == $3 {print $0}' file # outputs: cont1005 3642 3642 cont1037 3146 3146 cont1056 934 934 cont1059 1750 1750 cont1072 2577 2577 Note ...


6

I like to use the join command for this. For your example above, you can simply use the following: $ join file2 file1 ENST00000000412 0.779142 M6PR ENST00000001008 0.738143 FKBP4 ENST00000002501 0.715315 DBNDD1 ENST00000002596 0.713664 HS3ST1 This assumes your files are sorted. If they aren't sorted, you can use input redirection if your shell supports it (...


6

Invest some time in learning unix tools such as grep, sed and awk. sed -r '/^>/s/:[0-9]+-[0-9]+//' input.fasta > output.fasta This line of code matches the pattern :<numbers>-<numbers> on lines that start with > (FASTA sequence headers) and replaces that pattern with '', essentially removing text that matches the pattern.


6

Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match. zcat My_Hiseq_Data.fq.gz | \ grep -Eo '.{0,3}GATCGATC.*' | \ sed -En 's/.*/ \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \ grep --color=always GATCGATC Here are some ...


5

Assuming all your vcfs are in the same folder: for file in *.vcf; do grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt done


5

Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible. Here's a Perl script that should do exactly ...


5

awk -F "::" '{if($1~">"){gsub(">","");print ">"$2"_"$1} else {print $0}}' foo.fa Basically the same as from Devon but using -F to indicate the initial field delimiter and then using gsub to remove the >.


5

You have an obvious error in your shell script that should have given you an error message: sra_code = "DRR163""$i" Should be: sra_code="DRR163""$i" You cannot have spaces around the = sign in a variable assignment. The line sra_code = "DRR163""$i" means "run the command sra_code with the ...


5

This is an error from bwa saying that it can't find your input files: Here's a checklist: Check that you have no whitespace after your \ newline breaks. I.e. \ is the very last character Check that you have no whitespace in any path names Quote your variables: $REF_PATH -> "$REF_PATH" Use set -x at the top of your bash script to debug. In ...


4

multiqc kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the multiqc report and the actual fastqc html report to get a better idea of any error modes. Going off of @Kubator's answer, I noticed that there was no ...


4

The following sed command will do the trick. sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna Lines starting with > will be escaped and the lowercase characters will be replaced with - .


4

You should be able to parse out what you need using the tags in the .bam. 10xGenomics' website says what tags they add. https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -> fastq pipeline won't do ...


4

You could simply use bcftools sort for it: $ bcftools sort input.vcf > output.vcf If you really want to use bash only, you can do this: $ grep "^#" input.vcf > output.vcf $ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf The first command will write the header information to the new vcf file. The second will sort by contig name and ...


4

Go to the common parent folder under which all FASTQ.GZ files reside. They could be in sub-folders under the folder, but all of them should have the folder you're in as a common parent folder. From such a folder, you should be able to run a find find . -name "*.fastq.gz" to get your FASTQ.GZ files. Now, there are surely better ways of doing this but here's ...


4

The bug in your loop appears be "> filename" per column, this over-rides the last iteration resulting in a file with a single column, which has been over-written 85999 times. Thus the file output you observe should be the single column present in the last file of your array. If you replaced "> filename" with ">> filename" within the loop you ...


4

Reverse complement FASTA/Q: seqtk seq -r in.fa > out.fa https://github.com/lh3/seqtk


4

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP 'TTCG[CGT][ATGC]{1,15}TTCG[CGT]' [CGT] looks for one occurrence of one of the three. [ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases. Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors ...


4

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f. EDIT AS THE OP HAS EDITED THE QUESTION: First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the ...


4

First of all you are getting an error because you are missing a then in the third line of your code and then a fi to close the if. It should be: for i in *.vcf do if grep -q $i 1.txt; then cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/ fi done And I guess the one-liner below achieves what you are trying to do: while IFS= read -r line; do mv "$line"...


3

We sequence DNA on daily basis. Like terdon said, you can use cat to group all the FastQ of all your sample whether they are bgzipped or not like so: root@slurm1: ~ # ls am224_S6_L001_R1_001.fastq.gz am224_S6_L001_R2_002.fastq.gz am224_S6_L002_R1_003.fastq.gz am224_S6_L002_R2_004.fastq.gz am232_S10_L001_R1_001.fastq.gz am232_S10_L001_R2_002.fastq.gz ...


3

You could try it with sed. Replace all lower-case chars with - in lines not starting with >: sed -e '/^>/!s/[a-z]/-/g' in.fa


3

The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of sudo ability is not an issue. I have installed samtools and many other programs necessary for the workflow, however, once I am on the ...


3

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t"). In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter. egrep -v "/db_xref=\"taxon"| tr "\"" ...


3

You were so close! This should do the trick: for i in {314..727}; do sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna done Note that I am using -i which makes sed edit the original file. Because I use -i.bak, this will create a backup file with the original data. Once you are sure the command does what you need, you can delete all of these files ...


3

You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment. From bwa docs: If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair. BWA will not search through FASTQ files ...


3

You don't need a loop. You can do the whole thing with a simple awk one-liner: awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta I ran this on your files and got: $ for file in gene*; do echo "=== File: $file ==="; cat $file; done === File: gene1.fasta === >locus1 ATGCGTAGAG >ysr ATGTAGCGA >siv TAGTAGTAT === ...


Only top voted, non community-wiki answers of a minimum length are eligible