landsat - Image stacking with Python/Gdal

I am trying to accomplish some image stacking task in python. I am not proficient in Python and was hoping I can get some guidance to get started.

So I have large number of landsat images. Landsat images for one particular year are stored in a folder and there are sub-directories for each scene. For each scene there are numerous tif images for different VIs and I need to stack them according to their dates.I can stack images in one single folder by reading each image. But I want to automate the process as different path_rows have different number of images(one path_row might have 5 image but other one might have 10 image).

Answer

So this was my solution:

    test_vrt= ["gdalbuildvrt", "-separate", 'NDVI_test'+ ".vrt"]
    image_list =  glob.glob('*/*band1.tif')
      for myfile in image_list:
        test_vrt.append(myfile)
        subprocess.call(test_vrt)
    img2tiff =["gdal_translate","test.vrt", "test_composite.tif"]
    subprocess.call(img2tiff)