Hot answers tagged

22

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


9

If you want to avoid using extra libraries for any reason, you can just use a simple Python script (version 3.6 and above) to do this: fr = open("dup_test.fasta", "r") fw = open("dup_edited.fasta", "w") seq_dict = {} curr_header = '' for line in fr: line = line.strip() if line[0] == '>': if line not ...


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, ...


8

You can use seqkit for this purpose. seqkit rmdup -n seqs.fa -o seqs_without_duplicate.fa


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

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


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

I've been using this: cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c Explanation : paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read cut -f2 : print only the second column, to access the sequence after the paste wc -c : count the characters tr -d '\n': to remove from count 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

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"...


4

One way using awk: awk '/^>/ { f = !a[$0]++ } f' seqs.fa Explanation: The above is basically the same as: awk '/^>/ { f = !($0 in a); a[$0]++ } f' seqs.fa That is, for each header line (lines that start with a '>' character), set a flag, 'f' to true if the current line is not in an array, 'a'. Then add the (header) line to the array. Note that in ...


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

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

How to re-analyze 10X BAM files? This is a great question and honestly, I don't think there was an easy way to do this at the time the question was asked. The reason is that if you want to re-do the authors' analysis to get the gene expression (gene-cell / feature-barcode) matrix and they used cellranger (the official 10X pipeline software) to do that, you ...


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