How does the function $f(q)=\sum_{p\in\mathbb{P}} q^p$ behave near $q=1$? I wonder, is there an elementary or other closed-form function $g(q)$ which behaves like $$f(q)=\sum_{p\in\mathbb{P}} q^p$$ near $q=1$?
The plot of $f(q)$ looks like the following, but I failed to find an elementary function with a similar behavior near $q=1$:

 A: Here’s a perhaps slightly more pedestrian solution that leads to the concrete asymptotic behaviour of your function.
As $q\to1$, your sum averages over more and more primes, so asymptotically it doesn’t depend on the details of the distribution of the primes. Thus, we can replace the sum over the primes by an integral over their density $\frac1{\log p}$:
$$
\sum_pq^p\sim\int_a^\infty q^p\frac{\mathrm dp}{\log p}\;,
$$
where I left the lower limit $a$ unspecified; it doesn’t matter asymptotically and will drop out in due course. Now set $q=1-\epsilon$ to find
\begin{eqnarray}
\int_a^\infty q^p\frac{\mathrm dp}{\log p}
&=&
\int_a^\infty(1-\epsilon)^p\frac{\mathrm dp}{\log p}
\\
&=&
\int_a^\infty\mathrm e^{p\log(1-\epsilon)}\frac{\mathrm dp}{\log p}
\\
&\sim&
\int_a^\infty\mathrm e^{-\epsilon p}\frac{\mathrm dp}{\log p}
\\
&=&
\int_{\epsilon a}^\infty\mathrm e^{-x}\frac{\mathrm dx}{\epsilon\log\frac x\epsilon}
\\
&=&
\frac1\epsilon\int_{\epsilon a}^\infty\mathrm e^{-x}\frac{\mathrm dx}{\log x-\log\epsilon}
\\
&=&
-\frac1{\epsilon\log\epsilon}\left(\int_0^\infty\mathrm e^{-x}\mathrm dx+\frac1{\log\epsilon}\int_0^\infty\mathrm e^{-x}\log x\mathrm dx+O\left(\frac1{\log^2\epsilon}\right)\right)
\\
&=&
-\frac1{\epsilon\log\epsilon}\left(1-\frac{\gamma}{\log\epsilon}+O\left(\frac1{\log^2\epsilon}\right)\right)\;,
\end{eqnarray}
where $\gamma$ is the Euler–Mascheroni constant. Thus
$$
\sum_pq^p\sim-\frac1{(1-q)\log(1-q)}\;.
$$
To visualize the asymptotic behaviour, we can look at the reciprocal of your function:
\begin{eqnarray}
\left(\sum_pq^p\right)^{-1}
&=&
-\epsilon\log\epsilon\left(1-\frac{\gamma}{\log\epsilon}+O\left(\frac1{\log^2\epsilon}\right)\right)^{-1}\;,
\\
&=&
-\epsilon\left(\log\epsilon+\gamma+O\left(\frac1{\log\epsilon}\right)\right)^{-1}\;.
\end{eqnarray}
This shows why I carried along the term in $\epsilon\log^{-2}\epsilon$, even though it’s not part of the leading asymptotic behaviour, as it leads to a term linear in $\epsilon$ in the reciprocal that we need for a good comparison.
Here’s a graph of the reciprocal of your function (in purple) compared to $-(1-q)(\log(1-q)+\gamma)$ (in green):

A: It is part of the prime number theorem and partial summation that as $x\to 0,x >0$ $$\sum_p e^{-px} \sim \sum_{n\ge 2} \frac{e^{-nx}}{\log n}$$
The RH gives an error term $O(\sum_n n^{-1/2+\epsilon}e^{-nx})$.

The residue theorem explicit formula gives more information. 

As usual we look at the simpler series
$$\sum_n \Lambda(n)e^{-nx}=\frac1{2i\pi}\int_{2-i\infty}^{2+i\infty} \frac{-\zeta'(s)}{\zeta(s)}\Gamma(s)x^{-s}ds$$ $$=\sum Res(\frac{-\zeta'(s)}{\zeta(s)}\Gamma(s)x^{-s})=x^{-1}-\sum_\rho \Gamma(\rho) x^{-\rho}+f(x)+g(x)\log x$$ where $f,g$ are entire.
The latter series stays true by analytic continuation for $\Re(x) > 0$ using the density of zeros and the $\sqrt{\pi} e^{-\pi t/2}$ decay of $\Gamma(1/2+it)$ obtained from $\pi/\sin(\pi s)=\Gamma(s)\Gamma(1-s)$. When separating the $\Im(\rho)>0$ and $\Im(\rho)<0$ parts we get two functions that can be analytically continued along either clockwise or counterclockwise rotation around $0$. The same analysis works for $\sum_{p^k}\frac1k e^{-p^k x}$ and $\sum_p e^{-px}$ except that we'll get ugly $li(x^{-\rho})$ terms.

$f(x)+g(x)\log x$ corresponds more or less to $\sum_n n^{-1}\Lambda(n)e^{-x/n}$, this is the functional equation.

