Create PostGIS geometry tables from GTFS stops and shapes. Then join geom tables with route info. Tested with Corvallis GTFS data, which is the reason for reprojecting from EPSG:4326 to EPSG:2913. Also, SQL is weird.

gtfs_join_routes_to_stops_shape_geom.sql

create table route_lines as (			    
	select r.*, s.geom 
	from (		
			select route_id, route_short_name, route_long_name, route_color, min_headway_minutes
			from routes
			group by route_id
		) r
	join (
			select route_id, shape_id 
			from trips
			group by route_id, shape_id
		) t on t.route_id = r.route_id
	join 
		shape_lines s	on t.shape_id = s.shape_id
);

create table route_stop_pts as (
	select t.route_id, st.stop_id, s.geom
	from (
		select trip_id, route_id
		from trips
		group by trip_id, route_id
	) t
	join (
		select trip_id, stop_id
		from stop_times
		group by trip_id, stop_id
	) st on st.trip_id = t.trip_id
	join stop_pts s on s.stop_id = st.stop_id
	group by t.route_id, st.stop_id, s.geom
);

gtfs_stops_shapes_to_geom.sql

create table stop_pts as (
	select *, st_transform(st_setsrid(st_makepoint(stop_lon,stop_lat),4326),2913) as geom
	from stops
);
create index stop_pts_gist_idx on stop_pts using GIST(geom);
vacuum analyze stop_pts;

create table shape_lines as (
	with shape_pts as (
		select
			shape_id,
			shape_pt_sequence,
			st_transform(st_setsrid(st_makepoint(shape_pt_lon,shape_pt_lat),4326),2913) as geom
		from
			shapes
		order by
			shape_id, shape_pt_sequence
	)
	select
		shapes.shape_id,
		st_makeline(array_agg(shape_pts.geom order by shape_pts.shape_pt_sequence)) as geom
	from 
		shape_pts,
		shapes
	where
		shape_pts.shape_id = shapes.shape_id
	group by
		shapes.shape_id
);
create index shape_lines_gist_idx on shape_lines using GIST(geom);
vacuum analyze shape_lines;

results.md

Routes. Stops. Larger stop symbol = more routes at stop.