cannin/biostars_qa huggingface dataset:
Question: Hello Biostars,
Can anyone tell me how to prepare input data set for GSEA after Differential Gene Expression Analysis by DESeq2?
How will I rank the genes? Should I rank based on log2FC or Adjusted P value? Is there any way to generate a GSEA ready data directly from DESeq2?.
I was using topGo for gene ontology enrichment analysis before and recently came across GSEA. Which one is better GO enrichment analysis or GSEA? Even after going through the papers I couldn't find a significant difference between above two.

Thank you

Answer: Hi Sreeraj

Genes can be ranked based on fold change and P value and that can be used in GSEA package.

you can use this R code for this purpose.

x <- read.table("DE_genes.txt",sep = "\t",header = T)
head(x)
x$fcsign <- sign(x$log2.fold_change.)
x$logP=-log10(x$p_value)
x$metric= x$logP/x$fcsign
y<-x[,c("Gene", "metric")]
head(y)
write.table(y,file="DE_genes.rnk",quote=F,sep="\t",row.names=F)
Question: Hey biostars,

How I can align multiple fastq files on reference genome or I have to align them individually?

Answer: You can definitely align multiple FASTQ files to a reference genome - I do it all the time. You'd ideally be able to pass as parameters pairs of PE fastq files and any number of SE fastq files. Check out the documentation of the alignment program you're using to figure out the exact method to pass multiple files.
Question: Hi All,

**library(TCGAbiolinks)** is giving errors as *gridExtra* was renamed so I can not use it to download and analyzing the TCGA data. Except **GDC Portal** (our IT support having problems to install the data transfer tool, so I can not use it right now) is there any other alternatives to get the TCGA-BRCA data to do differential gene expression analysis.

Thanks,

Gokce

Answer: Can you show what you are doing and replicate the error so as to make up help you, if a dependency package or module is renamed then it should be also forked into the TCGAbiolonks package if all install.dependancies are made TRUE. If it's a bug then you can always contact their GitHub repo to report it.

Alternative is [TCGA2STAT][1]. Take a look at it. I believe you are just downloading the open files and not controlled access data so this package should also work and then you can use your custom lima-voom, edgeR or DESeq2 for your downstream analysis once you have all the meta-data and the count data downloaded for yourself.

[1]: https://cran.r-project.org/web/packages/TCGA2STAT/vignettes/TCGA2STAT.html
Question: Hello,

I was asked to find which CpGs are unique of EPIC 850K methylation array that are not observed by RRBS (reduced representation bisulfite sequencing). I have EPIC 850K manifest, but is there a public site where I can download all coordinates/locations for all CpGs detected by RRBS?

Do you think there will be CpGs that EPIC 850K covers that RRBS does not?

Thanks a lot for your help.

Answer: There won't be an exact list of CpGs covered by RRBS because it depends on the exact enzymes used and how tightly you perform size selection. I would propose that you perform the following procedure:

1. Use biopython to determine all possible fragments generated by the restriction enzymes you'll be using (there are some convenient functions for performing restriction digests on sequences in that package).
2. Determine a rough range of sequencable fragments, which will likely be something like 75-500 bases.
3. Choose a read length (N), because the results of all of this will be length-dependent.
4. For each of the fragments you selected from step 2, write the regions corresponding to the first/last N bases to a file in BED format.
5. Load the BED file from step 4 into an interval tree (there might be something in biopython for this, worst case scenario you can use deeptoolsintervals from deepTools).
6. Use biopython to iterate over the CpGs and query them for overlaps with the interval from step 5.
7. Write output files appropriately
8. Compare them to what the EPIC 850K covers.

Note that the EPIC 850K may give a ballpark estimate of all of this in their sales materials. I wouldn't be surprised if the EPIC 850K covers some CpGs that RRBS doesn't.
Question: I have gff files and files containing intergenic regions for different bacterial genomes. I want to extract the adjacent genes and their co-ordinate which corresponding those intergenic regions for all of these genomes. Are there any scripts or tools are available?

Example:

For *Bacillus anthracis*, I have two files:

1. gff file
2. files containing coordinates of intergenic regions:

`intergenic_region.txt`

```
185457-185562
320958-321064
1146951-1147049
1285399-1285500
3894344-3894451
4075706-4075815
```

I want to extract the coordinate of adjacent genes for each of those intergenic regions from gff files.

I have to do such kind of process for 300+ bacterial genomes.

Answer: Convert the GFF file to BED with BEDOPS *gff2bed*:

$ gff2bed < annotations.gff > annotations.bed

Convert the intergenic regions to BED. Assuming they are all on chromosome `chr1` and 0-indexed:

$ awk 'BEGIN { FS = "-"; } { print "chr1\t"$1"\t"$2; }' intergenic_region.txt > intergenic_region.bed

You may need to modify the name used for the chromosome, depending on the naming scheme used in the annotations file. You could browse through the first column of `annotations.bed` to see what the chromosome naming scheme looks like.

Once you have the annotations and intergenic regions in BED format, you can use BEDOPS *closest-features* to find the nearest annotations to the regions:

$ closest-features intergenic_region.bed annotations.bed > answer.bed
Question: I am using GAPIT package in R for GWAS. Already installed the package and import the genotype and phenotype file in R

myY <- read.table("mdp_traits.txt", head = TRUE), myG <- read.table("RabiGeno.hmp.txt", head = FALSE) and run the GAPIT using command
myGAPIT <- GAPIT(
Y=myY,
G=myG,
PCA.total=3

And the output is given below:

`There is an **Error in plot.window(...) : need finite 'xlim' values** .` So, how to solve this problem?

"--------------------- Welcome to GAPIT ----------------------------"
[1] "Converting genotype..."
[1] "Converting HapMap format to numerical under model of Middle"
[1] "Perform numericalization"
[1] "Succesfuly finished converting HapMap which has bits of 1"
[1] "Converting genotype done."
[1] "Calling prcomp..."
[1] "Creating PCA graphs..."
[1] "Joining taxa..."
[1] "Exporting PCs..."
[1] "PC created"
[1] "Calculating kinship..."
[1] "Number of individuals and SNPs are 182 and 45641"
[1] "Calculating kinship with VanRaden method..."
[1] "substracting P..."
[1] "Getting X'X..."
[1] "Adjusting..."
[1] "Calculating kinship with VanRaden method: done"
[1] "kinship calculated"
[1] "Creating heat map for kinship..."
[1] "Kinship heat map created"
[1] "Adding IDs to kinship..."
[1] "Writing kinship to file..."
[1] "Kinship save as file"
[1] "Kinship created!"
[1] "GAPIT.Genotype.View .Two pdf generate.successfully!"
[1] "--------------------Processing traits----------------------------------"
[1] "Phenotype provided!"
[1] "GAPIT.Phenotype.View in press..."
[1] "GAPIT.Phenotype.View .output pdf generate.successfully!"
[1] "GAPIT.Phenotype.PCA.View"

**Error in plot.window(...) : need finite 'xlim' values**

In addition: Warning messages:
1: In max(y2[, 2]) : no non-missing arguments to max; returning -Inf
2: In min(y2[, 2]) : no non-missing arguments to min; returning Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf
5: In min(x) : no non-missing arguments to min; returning Inf
6: In max(x) :
Show Traceback

Rerun with Debug
Error in plot.window(...) : need finite 'xlim' values

Answer: I've encountered this before - GAPIT needs numerical values for chromosome names to plot the Manhattan plot. If your chromosomes are named 'Chr1' or so it won't work, they have to be named 1, 2, 3, ....
Question: Can anyone share me a good GATK4 workflow, I Couldnt find anything good that explanins the process in a crisp and easy way.

Answer: If you're looking for a streamlined yet comprehensive explanation of GATK4 and the main Best Practices workflows, you may be interested in the O'Reilly book I co-authored called Genomics in the Cloud. (The cloud part is secondary, all the GATK-focused explanations are applicable pretty much regardless of where you will run the workflows).

See https://oreil.ly/genomics-cloud
Question: What is the meaning of **(-) log10(FDR)** value

How to interpret higher negative value to lower negative value by comparing each other?

e.g. 1:

Ribosome : *(-) log10(FDR) = 25*

Golgi : *(-) log10(FDR) = 5*

e.g. 2:
[![enter image description here][1]][1]

[1]: https://i.stack.imgur.com/zCAMa.jpg

Answer: It is just to make the value easier to read.

Instead of 1E-5 you are shown 5 and instead of 1E-25 you are shown 25

The bigger number means a smaller FDR
Question: I am working on Arabidopsis thaliana omics project. How can I get short and proper gene description for each AGI (Gene IDs) numbers?
Downloaded database from Tair10 has pretty lengthy names which are difficult to work with (plotting , summarizing in downs-stream works)

e.g.

AT1G01050	Soluble inorganic pyrophosphatase 1 OS=Arabidopsis thaliana (sp|q93v56|ipyr1_arath : 419.0)

AT1G01800	Enzyme classification.EC_1 oxidoreductases.EC_1.1 oxidoreductase acting on CH-OH group of donor(50.1.1 : 434.7) & (+)-neomenthol dehydrogenase OS=Arabidopsis thaliana (sp|q9m2e2|sdr1_arath : 357.0) (original description: none)

From where I can get / how to **`modify the names short and sweet`**?

Answer: - Go to ensembl's [BioMart][1]
- Choose `Dataset`: `Ensembl Plant Genes` and `Arabidopsis thaliana genes`
- Choose `Filters`->`Genes` -> `Input external references ID list` -> `Gene Stable ID(s)` and paste your IDs into the textfield
- Choose `Attributes`->`Gene` and select `Gene Stable ID` and ` Gene name`
- Click `Result` and download in the format you like

You can extract the IDs from the file example above with a simple `cut -f1 input_file > gene_ids.txt`.

fin swimmer

[1]: http://plants.ensembl.org/biomart/martview
Question: Hello everyone

I tried to install ggplots packages of R in Ubuntu server, but there were some errors which displayed as follows:

Warning messages:
1:In install.packages("ggplot2"):
Installation of package 'plyr' had non-zero exit status

2: In install.packages("ggplot2"):
Installation of package 'reshape2' had non-zero exit status

3: In install.packages("ggplot2"):
Installation of package 'scales' had non-zero exit status

4: In install.packages("ggplot2"):
Installation of package 'ggplot2' had non-zero exit status

Then I tried to use ggplot2 through command:

library('ggplot2')

But terminal told me that

Error in library (ggplot2): there is no package called 'ggplot2'

I have found many information about this problem on internet, but no available ways find and error was error.

If you have some suggestion, please tell me.

Thanks for your help.

Answer: It seems that you are using a very old version R; so, even if you succeed in installing ggplot2, be aware that it may not work correctly afterwards. In any case, to solve the issue you should install Rcpp separately - e.g. `sudo apt-get install r-cran-rcpp`
Question: Hey all :]

I use samtools's depth, and occasionally samtool's pileup commands to calculate coverage of my reads to the genome before binning for coverage. I'm pretty sure everyone else does too ;)

One very common problem is that people find "samtools depth" and "samtools mpileup" don't match up - and it's commonly attributed to filtering poor quality reads, duplicates, etc (mpileup doing more filtering).

But there's a ton of other questions I just can't get the answers too from the sametools docs, namely:

Does depth/pileup count the region between a paired reads, or just the reads themselves?

Does samtools depth/pileup count singletons in paired end sequencing? Does it extrapolate based on average read length to fillout the whole fragment?

If a read maps to multiple locations, is it counted multiple times?

I'm working with ChIP-Seq data, so i might have to correct for peak-shift. What do you guys think?

Thanks! :)

Answer: There comes a time in every bioinformatician's life when he/she needs to just read the source code. This is one of those times.

So, depth actually uses pileup internally, just to confuse things a bit more. Depth does not count regions spliced over (i.e., covered by an N operator) or between reads in a pair, which is good because what's actually going on there would be undefined. I don't think that pileup counts regions between reads, though I guess I'd have to look to be 100% certain.

Neither tool cares whether a read is a singleton or a part of a pair unless you tell them to filter according to that. Samtools is a pretty general program, so it won't do ChIP-seq specific things like extrapolating bounds based on average template length.

In general, you'll need to start going through the code to get answers to questions like these.
Question: This is cross-posted [from here](https://bioinformatics.stackexchange.com/q/4419/351).

I would like to use uniprot's [sparql endpoint](https://sparql.uniprot.org/) to retrieve all proteins that

1. are reviewed (required)
2. are associated with taxonomy IDs `562` and `3702` (required)
3. have a KO associated with them (optional)
4. "evidence for the existence of a protein " should be either on protein or transcript level (required)
5. have an EC number associated with them (required)

I have so far (points 1 and 2):

PREFIX up:
PREFIX taxon:
PREFIX rdfs:

SELECT ?protein ?taxon ?name
WHERE
{
?taxon a up:Taxon .
?taxon up:scientificName ?name .
VALUES ?taxonlist { taxon:562 taxon:3702 }
?taxon rdfs:subClassOf ?taxonlist .

?protein a up:Protein .
?protein up:organism ?taxon .
?protein up:reviewed true . # have to be reviewed

}

This, however, does not return anything for `3702`. How can this be fixed and how can I incorporate points 3-5?

Additionally, is there now a way to connect uniprot's sparql endpoint with rhea's [sparql endpoint](https://sparql.rhea-db.org/sparql) to retrieve all associated reactions and their stoichiometries (with ChEBI IDs) for the selected proteins from above? Example 19 seems to suggest that this connection is possible but I am not quite sure how to accomplish it.

Answer: 1) Is correct in the query with

?protein up:reviewed true .

2) The query in the Q. does not return anything for taxon:3702 as there are no rdfs:subClasses for Aribidopsis Thaliana, it is a leaf node. This means the entry is directly linked to that taxon instead of via it's ancestors. This is fixed by changing the query slightly to deal with both the ancestor and direct case (both sides of the UNION below)

VALUES ?taxonlist { taxon:3702 taxon:562}
{
?taxon rdfs:subClassOf ?taxonlist .
?protein up:organism ?taxon .
} UNION {
?protein up:organism ?taxonlist .
}

3) we use the cross reference section which are done via rdfs:seeAlso . But as there is the possibility of more than one KO per entry we group them with a subquery.

OPTIONAL {
SELECT ?protein (GROUP_CONCAT(?ko; SEPARATOR=", ") AS ?kos)
WHERE{
?protein rdfs:seeAlso ?ko .
?ko up:database
PREFIX taxon:
PREFIX rdfs:

SELECT
?protein
?taxon
?name
?kos
?ecs
WHERE
{
?protein a up:Protein .
?protein up:reviewed true . # have to be reviewed
?taxon a up:Taxon .
?taxon up:scientificName ?name .
VALUES ?taxonlist { taxon:3702 taxon:562 }
{
?taxon rdfs:subClassOf ?taxonlist .
?protein up:organism ?taxon .
} UNION {
?protein up:organism ?taxonlist .
}

{
?protein up:existence up:Evidence_at_Protein_Level_Existence .
} UNION {
?protein up:existence up:Evidence_at_Transcript_Level_Existence .
}
{
SELECT ?protein (GROUP_CONCAT(?ec; SEPARATOR=", ") AS ?ecs)
WHERE{
?protein up:enzyme|((up:component|up:domain)/up:enzyme) ?ec
} GROUP BY ?protein
}
OPTIONAL {
SELECT ?protein (GROUP_CONCAT(?ko; SEPARATOR=", ") AS ?kos)
WHERE{
?protein rdfs:seeAlso ?ko .
?ko up:database
} GROUP BY ?protein
}
}

Which is testable at s[parql.uniprot.org][1].

[1]: http://tinyurl.com/y7mv68gp
Question: Terminology suggestions?

Original question removed, since it does not make sense. Changing it so that Kao's answer fits.

Answer: I've heard of the term "epiallels".

http://www.nature.com/pr/journal/v61/n5-2/full/pr2007128a.html
Question: Dear all,

I am new to the field. I am trying to analyze single end 100b FastQ files with ~70million reads/sample. I am trying to determine if adapter sequences are present and if so how to go about them. I ran FastQC on the files and reports show they each have an "overrepresented sequence" of an "illumina index adapter" in them.

I have the following questions:

1. Does sample1 look like a trimmed file or it requires adapter trimming?

2. If further trimming is recommended what would be the best seq/adapter option to be used for cutadapt/TrimGalore? [See below for my thoughts so far]

3. Based on the FastQC report, do I need to worry about presence of any other adapter sequences beside the index?

My thoughts on question 2:
The sequences for illumina index adapter format appear to be:

GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG

These are the adapter sequences found in my FastQC report for sample 1:

GATCGGAAGAGCACACGTCTGAACTCCAGTCACCATGGCATCTCGTATGC
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACCATGGCATCTCGTATG

I am thinking of using below options for cutadapt/trimgalore to remove the adapter(s):

trim_galore sample1.fastq.gz -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -q 20 --length 20 –fastqc

However, it seems that trimmomatics for instance only takes care of the initial sequence of the index adapter (only up to Ns and not after):
https://github.com/timflutre/trimmomatic/blob/master/adapters/TruSeq3-SE.fa

Many thanks for your time and reply beforehand.

Answer: 0.5%? You really don't have to worry about that if you don't want to.

The N's are for the variable index region. You know what the index is, you can see it in the fastqc report. Why would you put N's in?
Question: Hi!

Does anyone know any database, which collects a list of parasitic worm species? I've only found some dbs which describe diseases, caused by these organisms, and could not find a referenceable, well-organized list of them anywhere.

Thanks!

Answer: there is a **Species interaction datasets** database available for download at https://www.globalbioticinteractions.org/data.html

$ wget -O - -q "https://depot.globalbioticinteractions.org/snapshot/target/data/tsv/interactions.tsv.gz" | gunzip -c | cut -f1,3,33,35,37 | grep -i parasite

EOL:12001247	Leptoconchus massini	parasiteOf	EOL:200339	Fungia repanda
EOL:12001247	Leptoconchus massini	parasiteOf	EOL:200324	Fungia concinna
EOL:12001243	Leptoconchus inpleuractis	parasiteOf	EOL:200337	Fungia paumotensis
EOL:12001243	Leptoconchus inpleuractis	parasiteOf	EOL:200336	Fungia moluccensis
EOL:12001243	Leptoconchus inpleuractis	parasiteOf	EOL:200333	Fungia gravis
EOL:12001242	Leptoconchus inpileus	parasiteOf	EOL_V2:200352	Halomitra pileus
EOL:12001241	Leptoconchus inlimax	parasiteOf	EOL_V2:200354	Herpolitha limax
EOL:12001239	Leptoconchus ingranulosa	parasiteOf	EOL:200332	Fungia granulosa
EOL:12001238	Leptoconchus ingrandifungi	parasiteOf	EOL_V2:200365	Sandalolitha dentata
EOL:12001237	Leptoconchus infungites	parasiteOf	FBC:SLB:SpecCode:45838	Fungia fungites
EOL:12001236	Leptoconchus incycloseris	parasiteOf	EOL:200325	Fungia costulata
WORMS:847839	Umimayanthus parasiticus	parasiteOf	WORMS:166232	Callyspongia (Cladochalina) vaginalis
WORMS:847839	Umimayanthus parasiticus	parasiteOf	WORMS:166232	Callyspongia (Cladochalina) vaginalis
EOL_V2:620727	Spirobranchus giganteus	parasiteOf	EOL_V2:538613	Neofibularia nolitangere
EOL_V2:620727	Spirobranchus giganteus	parasiteOf	EOL_V2:335778	Ircinia felix
EOL_V2:620727	Spirobranchus giganteus	parasiteOf	WORMS:169016	Monanchora arbuscula
EOL_V2:620727	Spirobranchus giganteus	parasiteOf	WORMS:169121	Desmapsamma anchorata
EOL_V2:620727	Spirobranchus giganteus	parasiteOf	WORMS:166232	Callyspongia (Cladochalina) vaginalis
WORMS:328023	Fauchaldius cyrtauloni	parasiteOf	WORMS:171796	Cyrtaulon sigsbeei
WORMS:421835	Upogebia annae	parasiteOf	WORMS:131849	Xestospongia sp.
EOL_V2:1039933	Leucothoe spinicarpa	parasiteOf	EOL_V2:335032	Spheciospongia vesparium

you'll 'just' need to determine if an organism is a worm or not.
Question: I have 30 small fastq files from same sample, and I want to merge it into one file. I know the command is

cat file1.fastq file2.fastq > bigfile.fastq

but is there any short cut for doing it? It just looks silly to type 30 file names one by one...

Thank you for any idea!

Answer: cat file*.fastq > bigfile.fastq
Question: I have a table:

W X Y Z
A 2 3 4
A 2 3 6
B 1 2 3
C 3 2 1
B 1 3 4
B 1 2 2

# reproducible data
table <- read.table(textConnection("W\tX\tY\tZ
A\t2\t3\t4
A\t2\t3\t6
B\t1\t2\t3
C\t3\t2\t1
B\t1\t3\t4
B\t1\t2\t2"), header = TRUE)

I want to combine all rows with same value in Column W, so that it looks like

W X Y Z
A 2 3 2,6
B 1 2,3 3,4,2
C 3 2 1

**Note:** All different values of column **W** entries appear as comma separated.

Answer: You're looking for the R package 'plyr'

X1 <- read.table(textConnection("W\tX\tY\tZ
A\t2\t3\t4
A\t2\t3\t6
B\t1\t2\t3
C\t3\t2\t1
B\t1\t3\t4
B\t1\t2\t2"),header=T)

library(plyr)

X2 <- ddply(X1, .(W), summarize,
Xc=paste(X,collapse=","),
Zc=paste(Z,collapse=",") ,
Yc= paste(Y,collapse=","))

X2

W Xc Zc Yc
1 A 2,2 4,6 3,3
2 B 1,1,1 3,4,2 2,3,2
3 C 3 1 2

X2u <- ddply(X1, .(W), summarize,
Xc=paste(unique(X),collapse=","),
Zc=paste(unique(Z),collapse=",") ,
Yc= paste(unique(Y),collapse=","))

X2u
W Xc Zc Yc
1 A 2 4,6 3
2 B 1 3,4,2 2,3
3 C 3 1 2
Question: I have a `VCF` file with `SNPs` and genes subset of `GFF` file (only genes are present). How to extract `SNPs` in `VCF` format located in genes from my data?

Answer: Use [bedtools][1]:

$ bedtools intersect -a input.vcf -b genes.gff -header -wa > output.vcf

---
**EDIT:**

For (very) large `vcf` files it might be more efficient to `bgzip` and `tabix` index the `vcf` file, convert your `gff` to `bed` and use `tabix` to query the regions

**1. bgzip and index**

$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

**2. gff to bed**

E.g with [BEDOPS][2]:

$ gff2bed < genes.gff > genes.bed

**3. Query the regions**

$ tabix -R genes.bed -h input.vcf.gz > output.vcf

fin swimmer

[1]: https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html?highlight=intersect
[2]: https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/gff2bed.html
Question: I've been using this method to rename files at the end of a step:

processed_fastq_1:
type: File
outputBinding:
glob: ${ return '**/' + inputs.fastq1.basename }
outputEval: |
${
self[0].basename = inputs.add_rg_SM + '_R1.fastq.gz';
return self[0]
}

But I'm not sure if this is an acceptable method. Should there be an intermediate step using python or calling to a script to rename files in the middle of a CWL workflow?

Answer: I would suggest adding a simple `mv` command to the end of the command line, via arguments, e.g.

```
arguments:
- position: 1000
shellQuote: false
valueFrom: ' && mv $(inputs.fastq1.basename) $(inputs.add_rg_SM)_R1.fastq.gz'
```
then you could just glob `$(inputs.add_rg_SM)_R1.fastq.gz`

This way you don't create duplicate files, which would be the case when using a separate step. Also, other users using the tool can see that renaming is happening just by looking at the command-line.
Question: I have two files.One file contains blast result of contigs and the next file contains the contigs that are used for blast.I want an output text file of sequences in a separate folder which gives no hits in blast result from original fasta file with all contigs. I have no idea about this.Can you please help me to create perl script for this problem.

Answer: 1) Extract names of contigs without hit in blast:

grep -B5 "***** No hits found" blast.txt | grep Query | sed 's/Query= //g' > ids

2) Extract fasta records of those contigs without hit:

cat ids | xargs -n 1 samtools faidx contigs.fa > contigs_nohit.fa