What are the modular forms which "naturally" corresponds to Riemann Zeta function or Dirichlet L function ? We know that, if we take the coefficients of the Dirichlet series of Riemann zeta function or Dirichlet L functions and use those coefficients in q-expansion, we will not get modular forms.
But is there any other ways which can "naturally" generate modular forms from 
Riemann zeta function or Dirichlet L functions ?
 A: There is a simple reason why you will never get the coefficients of the Dirichlet series $\zeta(s)$ or $L(s, \chi)$ as the q-expansion coefficients of a modular form. The reason is that these $L$-series have Euler products where the term for a prime $p$ is $1 / $(linear polynomial in $p^{-s}$), while modular forms always give you $1 / $(quadratic polynomial in $p^{-s}$). 
Of course, there's two easy ways to make quadratic polynomials out of linear polynomials:


*

*You can multiply two linear polynomials together. This leads you to the conclusion that products like $\zeta(s)^2$ or $L(s, \chi_1) L(s, \chi_2)$ should correspond to modular forms; and this is indeed the case, at least half of the time (there is a parity condition you have to impose). The corresponding modular forms are called Eisenstein series. An example is given in user1952009's answer, which constructs a modular form corresponding to the product $\zeta(s) \zeta(s + 1 - 2k)$. (The fact that $1-2k$ is odd is an instance of the parity condition: $\zeta(s) \zeta(s + 2)$ doesn't correspond to a modular form.) 

*You can square the variable in your linear polynomial to get a quadratic polynomial. This leads you to the conclusion that e.g. $\zeta(2s)$ should correspond to a modular form (because $1 - p^{-2s}$ is a quadratic polynomial in $p^{-s}$). The corresponding modular forms are called $\theta$-series. In fact, Riemann used this in his original 1859 paper where he introduced the Riemann hypothesis: his proof of the functional equation of the Riemann zeta function relies on the fact (due to Jacobi) that the $\theta$-series $\tfrac{1}{2} + \sum_{n \ge 1} q^{n^2}$ is a modular form of weight 1/2. 

A: $\scriptstyle \color{grey}{\text{(eliding the convergence problems)}}$
The Eiseinstein series $$G_{2k}(\tau) = \sum_{(n,m) \in \mathbb{Z}^2 \setminus (0,0)} \frac{1}{(m+n\tau)^{2k}} \tag{1}$$ Has the Fourier series representation $$G_{2k}(\tau) = 2 \zeta(2k) + 2 \frac{(2i\pi)^{2k}}{(2k-1)!} \sum_{n=1}^\infty \sigma_{2k-1}(n) e^{2i \pi n \tau}$$
Using $(2\pi n)^{-s} \Gamma(s) = \int_0^\infty x^{s-1} e^{-2 \pi nx}dx$
we have that the Mellin transform of $G_{2k}(ix)-2\zeta(2k)$ is 
$$F(s)  = \int_0^\infty x^{s-1} (G_{2k}(ix)-2\zeta(k))dx = 2\frac{(2i\pi)^{2k}}{(2k-1)!} \int_0^\infty x^{s-1}  \sum_{n=1}^\infty \sigma_{2k-1}(n) e^{-2 \pi n x}dx$$
$$ = 2\frac{(2i\pi)^{2k}}{(2k-1)!}\sum_{n=1}^\infty \sigma_{2k-1}(n)\int_0^\infty x^{s-1}e^{-2 \pi n x}dx = 2\frac{(2i\pi)^{2k}}{(2k-1)!} \Gamma(s) (2\pi)^{-s}\sum_{n =1}^\infty \sigma_{2k-1}(n) n^{-s}$$
$$ =  \boxed{\frac{2(-1)^{k}}{(2k-1)!} (2i\pi)^{2k-s}\Gamma(s) \zeta(s) \zeta(s-2k+1)}$$

Since $\zeta(s) = \sum_{n=1}^\infty n^{-s}$ and $\sigma_{2k-1}(n) = \sum_{d | n} d^{2k-1}$ we have $$\zeta(s) \zeta(s-2k+1) = (\sum_{l=1}^\infty l^{-s})(\sum_{m=1}^\infty m^{2k-1} m^{-s})=  \sum_{n=1}^\infty n^{-s} \sum_{d | n} d^{2k-1} =   \sum_{n=1}^\infty n^{-s} \sigma_{2k-1}(n)$$

Using $i \pi - 2 i \pi \sum_{d=0}^\infty e^{2 i \pi d \tau} =  i\pi - 2 i \pi \frac{1}{1-e^{2i \pi \tau}} = i \pi \frac{e^{2i \pi \tau}+1}{e^{2i \pi \tau}-1}=\pi \cot \pi \tau = \frac{1}{\tau} + \sum_{m=1}^\infty \frac{1}{\tau-m}+\frac{1}{\tau-m}$ and differentiating $2k-1$ times we have $$\frac{(- 2 i \pi)^{2k}}{(2k-1)!} \sum_{d=1}^\infty d^{2k-1} e^{2 i \pi d \tau} = \sum_{m=-\infty}^\infty  \frac{1}{(m+\tau)^{2k}}$$
So  that
$$G_{2k}(\tau) = \sum_{m \in \mathbb{Z}^*} \frac{1}{m^{2k}} + \sum_{n=1}^\infty \sum_{m=-\infty}^\infty \frac{1}{(m+n\tau)^{2k}}$$
 $$ = 2\zeta(2k) + \sum_{n=1}^\infty \frac{(- 2 i \pi)^{2k}}{(2k-1)!} \sum_{d=1}^\infty d^{2k-1} e^{2 i \pi d n \tau}$$
$$ = 2\zeta(2k)+\frac{(- 2 i \pi)^{2k}}{(2k-1)!}  \sum_{n=1}^\infty \sigma_{2k-1}(n) e^{2 i \pi  n \tau}$$

Of course, the functional equation for $\zeta(s)$ proves that $G_{2k}(\tau)$ is a modular form, but $(1)$ proves directly that $$G_{2k}(-1/\tau) = \tau^{2k} G_{2k}(\tau)$$
i.e. $G_{2k}(i/x)-2\zeta(2k) = (ix)^{2k} G_{2k}(ix) -2\zeta(2k)$
and
$$F(s)  = \int_0^\infty x^{s-1} (G_{2k}(ix)-2\zeta(2k))dx \underset{y = 1/x}{=} \int_0^\infty y^{-s-1} (G_{2k}(i/y)-2\zeta(2k))dy $$ 
$$ = \int_0^\infty y^{-s-1} ((iy)^{2k} G_{2k}(iy) -2\zeta(2k))dy = (-1)^{k} F(2k-s)$$
