I think you could be looking for expression quantitative trait loci (eQTL).
Here is a short overview article.
But I am not sure if that concept is applicable to your condensed data. You would need the original genotyping information to calculate an association test between variants and gene expression phenotype, not just whether there is any "interesting" (by what means?) variant in a gene (e.g. a SNP). Instead, you have to let the association test decide about variants that may be linked to gene expression. So the answer to Q1 is "yes, association testing".
There is an R-package called Matrix eQTL you could try. If you can convert your data into its input format, that should be fine, but I don't know if you have enough samples to achieve good power to detect anything.
For Q2 and Q3 you needed to define "background expression" first. In principle you could use your binary variable "hit/no hit" to divide samples into two groups (per gene) and do a test for differential expression. What is best for that purpose depends on how the expression values were obtained. A wilcoxon.test would always work but may be severely underpowered with your example data. On the other hand, standard DGE methods like DESeq2 use sample groupings that are the same for each gene and are therefore not applicable.
Edit: here is some R-code, thanks to the reproducible example you gave, that calculates t.test and wilcox.test for each gene, based on the grouping by phenotype.
set.seed(1234)
exp_matrix <- replicate(10, rnorm(4))
rownames(exp_matrix) <- paste0("gene_", letters[1:dim(exp_matrix)[1]])
colnames(exp_matrix) <- paste0("sample_", letters[1:dim(exp_matrix)[2]])
phe_martix <- replicate(10, ifelse(sample(16, 4) > 2, 0, 1))
rownames(phe_martix) <- paste0("gene_", letters[1:dim(exp_matrix)[1]])
colnames(phe_martix) <- paste0("sample_", letters[1:dim(exp_matrix)[2]])
# could also filter the matrix to only include rows with > 2 integration events
mat <- cbind(exp_matrix, phe_martix)
# run the combined matrix through the tests, use try because of lacking samples
ret <- apply(mat, 1, function(x){
exp <- x[1:ncol(exp_matrix)]
phe <- as.factor(x[(ncol(exp_matrix)+1):(2*ncol(exp_matrix))])
list (
t.test=try(t.test(exp ~ phe)),
w.test=try(wilcox.test(exp ~ phe)))
})
## returns a list of the test results for each gene
